C     Last change:  ERB  22 Aug 2000   10:31 am
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     AMG1R5                                        MAIN SUBROUTINE
C
C     RELEASE 1.5, OCTOBER 1990
C
C     CHANGES AGAINST VERSION 1.1, JULY 1985:
C
C 1.  A BUG WAS DETECTED WHICH UNDER CERTAIN CIRCUMSTANCES INFLUENCED
C     SLIGHTLY THE CONVERGENCE RATE OF AMG1R1. FOR THAT REASON, THE
C     FOLLOWING LINE IN SUBROUTINE RESC:
C     IW(IMAXW(KC-1)+1) = IA(IMIN(KC))
C     HAS BEEN CHANGED TO:
C     IW(IMAXW(KC-1)+1) = IAUX
C
C 2.  A BUG WAS DETECTED IN SUBROUTINE PWINT. UNDER CERTAIN CIRCUM-
C     STANCES AN UNDEFINED VARIABLE WAS USED. ALTHOUGH THIS DID NOT
C     AFFECT THE NUMERICAL RESULTS, PROBLEMS CAN OCCUR IF CHECKING
C     FOR UNDEFINED VARIABLES IS USED. TO FIX THIS ERROR, IN PWINT
C     THE LABEL 1000 WAS MOVED TO THE STATEMENT
C     IBLCK1 = IMINW(K).
C
C 3.  A PARAMETER LRATIO HAS BEEN INTRODUCED, DENOTING THE RATIO
C     OF SPACE OCCUPIED BY A DOUBLE PRECISION REAL VARIABLE AND
C     THAT OF AN INTEGER. FOR THE IBM-VERSION LRATIO HAS BEEN SET
C     TO 2. CHANGE THIS VALUE IF NECESSARY. (IF, FOR EXAMPLE, YOU
C     WANT TO CHANGE THE DOUBLE PRECISION VECTORS TO SINGLE PRE-
C     CISION, LRATIO HAS TO BE SET TO 1. IN THE YALE SMP - ROUTINE
C     NDRV THERE IS A PARAMETER LRATIO, TOO.
C
C 4.  TYPE DECLARATIONS REAL*4 AND REAL*8 HAVE BEEN CHANGED TO THE
C     STANDARD-CONFORMING KEYWORDS REAL AND DOUBLE PRECISION, RESPEC-
C     TIVELY.
C
C 5.  CALLS TO THE FOLLOWING INTRINSIC FUNCTIONS HAVE BEEN REPLACED BY
C     CALLS USING GENERIC NAMES: DSQRT, MIN0, MAX0, IABS, DABS, FLOAT,
C     DFLOAT, DMAX1, ISIGN, IDINT, DLOG10.
C
C 6.  A SAVE STATEMENT HAS BEEN INSERTED IN ALL SUBROUTINES.
C
C 7.  EXTERNAL DECLARATION STATEMENTS HAVE BEEN INSERTED IN ALL SUB-
C     ROUTINES FOR ALL EXTERNAL REFERENCES.
C
C-----------------------------------------------------------------------
C
C     CHANGE AGAINST VERSION 1.3, APRIL 1986:
C
C 1.  A BUG IN SUBROUTINE CHECK HAS BEEN REMOVED. IF THE ORIGINAL MATRIX
C     WAS STORED IN AN UNSYMMETRIC WAY, THE SYMMETRIZATION BY AMG1R3
C     COULD FAIL UNDER CERTAIN CIRCUMSTANCES. FOR A FIX, THE FOLLOWING
C     STATEMENTS IN SUBROUTINE CHECK HAVE BEEN CHANGED:
C
C     DO 450 J=IA(I)+1,IA(I+1)-1 WAS CHANGED TO
C     DO 450 J=IA(I)+1,ICG(I)-1
C
C     DO 430 J1=IA(I1)+1,IA(I1+1)-1 WAS CHANGED TO
C     DO 430 J1=IA(I1)+1,ICG(I1)-1
C
C     DO 550 J=IA(I)+1,IA(I+1)-1 WAS CHANGED TO
C     DO 550 J=IA(I)+1,ICG(I)-1
C
C     DO 530 J1=IA(I1)+1,IA(I1+1)-1 WAS CHANGED TO
C     DO 530 J1=IA(I1)+1,ICG(I1)-1
C
C 2.  THE EXPLANATORY PART IN SUBROUTINE AMG1R5 HAS BEEN ENLARGED TO
C     AVOID MISUNDERSTANDINGS IN THE DEFINITION OF THE ARGUMENT LIST.
C
C-----------------------------------------------------------------------
C
C     CHANGE AGAINST VERSION 1.4, OCTOBER, 1990 (BY JOHN W. RUGE)
C
C 1.  A BUG IN SUBROUTINE CHECK HAS BEEN REMOVED. IF THE ORIGINAL MATRIX
C     WAS STORED IN AN UNSYMMETRIC WAY, THE SYMMETRIZATION BY AMG1R3
C     COULD STILL FAIL UNDER CERTAIN CIRCUMSTANCES, AND WAS NOT FIXED
C     IN THE PREVIOUS VERSION. IN ADDITION, THE ROUTINE WAS CHANGED
C     IN ORDER TO AVOID SOME UNNECESSARY ROW SEARCHES FOR TRANSOSE
C     ENTRIES.
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
      SUBROUTINE AMG1R5(A,IA,JA,U,F,IG,
     +                  NDA,NDIA,NDJA,NDU,NDF,NDIG,NNU,MATRIX,
     +                  ISWTCH,IOUT,IPRINT,
     +                  LEVELX,IFIRST,NCYC,EPS,MADAPT,NRD,NSOLCO,NRU,
     +                  ECG1,ECG2,EWT2,NWT,NTR,
     +                  IERR)
C
C         -----------------------------------------------
C         | AMG-MODULE FOR SOLVING LINEAR SYSTEMS L*U=F |
C         -----------------------------------------------
C
C         --------------------------------------------------------------
C
C     ASSUMPTIONS ON L:
C
C         THE PROGRAM REQUIRES:
C
C             - DIAGONAL ENTRIES ARE ALWAYS POSITIVE (ON ALL GRIDS);
C             - L IS A SQUARE MATRIX WHICH IS EITHER REGULAR OR SINGULAR
C               WITH ROWSUMS=0.
C
C         FOR THEORETICAL REASONS THE FOLLOWING SHOULD HOLD:
C
C             - L POSITIVE DEFINITE (OR SEMI-DEFINITE WITH ROWSUM=0)
C             - L "ESSENTIALLY" POSITIVE TYPE, I.E.,
C
C                  -- DIAGONAL ENTRIES MUST BE > 0 ;
C                  -- MOST OF THE OFF-DIAGONAL ENTRIES <= 0 ;
C                  -- ROWSUMS SHOULD BE >= 0 .
C
C
C     THE USER HAS TO PROVIDE THE MATRIX L, THE RIGHT HAND SIDE F AND
C     CERTAIN POINTER VECTORS IA AND JA.
C
C         --------------------------------------------------------------
C
C     STORAGE OF L:
C
C         THE NON-ZERO ENTRIES OF THE MATRIX L ARE STORED IN
C         "COMPRESSED" SKY-LINE FASHION IN A 1-D VECTOR A, I.E., ROW
C         AFTER ROW, EACH ROW STARTING WITH ITS DIAGONAL ELEMENT. THE
C         OTHER NON-ZERO ROW ENTRIES FOLLOW THEIR DIAGONAL ENTRY IN ANY
C         ORDER.
C
C
C         IN ORDER TO IDENTIFY EACH ELEMENT IN A, THE USER HAS TO
C         PROVIDE TWO POINTER ARRAYS IA AND JA. IF NNU DENOTES THE TOTAL
C         NUMBER OF UNKNOWNS, THE NON-ZERO ENTRIES OF ANY ROW I OF L
C         (1.LE.I.LE.NNU) ARE STORED IN A(J) WHERE THE RANGE OF J
C         IS GIVEN BY
C
C                     IA(I) .LE. J .LE. IA(I+1)-1.
C
C         THUS, IA(I) POINTS TO THE POSITION OF THE DIAGONAL ENTRY OF
C         ROW I WITHIN THE VECTOR A. IN PARTICULAR,
C
C                     IA(1) = 1 ,  IA(NNU+1) = 1 + NNA
C
C         WHERE NNA DENOTES THE TOTAL NUMBER OF MATRIX ENTRIES STORED.
C         THE POINTER VECTOR JA HAS TO BE DEFINED SUCH THAT
C         ANY ENTRY A(J) CORRESPONDS TO THE UNKNOWN U(JA(J)), I.E.,
C         JA(J) POINTS TO THE COLUMN INDEX OF A(J).
C         IN PARTICULAR, A(IA(I)) IS THE DIAGONAL ENTRY OF ROW I
C         AND CORRESPONDS TO THE UNKNOWN U(I): JA(IA(I))=I.
C
C         IN THIS TERMINOLOGY, THE I-TH EQUATION READS AS FOLLOWS
C         (FOR ANY I WITH  1.LE.I.LE.NNU):
C
C                  F(I) =        SUM      A(J) * U(JA(J))
C                           J1.LE.J.LE.J2
C
C         WHERE F(I) DENOTES THE I-TH COMPONENT OF THE RIGHT HAND
C         SIDE AND
C
C                     J1 = IA(I) ,  J2 = IA(I+1)-1.
C
C         NOTES: THE ENTRY IA(NNU+1) HAS TO POINT TO THE FIRST FREE
C                ENTRY IN VECTORS A AND JA, RESPECTIVELY. OTHERWISE,
C                AMG CANNOT KNOW THE LENGTH OF THE LAST MATRIX ROW.
C
C                THE INPUT VECTORS A, IA AND JA ARE CHANGED BY AMG1R5.
C                SO, AFTER RETURN FROM AMG1R5, THE PACKAGE MUST NOT
C                BE CALLED A SECOND TIME WITHOUT HAVING NEWLY DEFINED
C                THE INPUT VECTORS AND USING ISWTCH=4. OTHERWISE, THE
C                SETUP PHASE WILL FAIL.
C                  ON THE OTHER HAND, RUNNING AMG A SECOND TIME ON THE
C                SAME INPUT DATA WITH ISWTCH=4 HAS NO SENSE, BECAUSE
C                THE RESULTS OF THE FIRST SETUP PHASE ARE STILL STORED
C                AND THUS THIS PHASE CAN BE SKIPPED IN A SECOND CALL.
C                IN ORDER TO DO THIS, SET ISWTCH TO 1, 2 OR 3.
C
C-----------------------------------------------------------------------
C
C         THE FORM OF THE CALLING PROGRAM HAS TO BE AS FOLLOWS:
C
C               PROGRAM DRIVER
C         C
C               DOUBLE PRECISION A(#NDA),U(#NDU),F(#NDF)
C               INTEGER IA(#NDIA),JA(#NDJA),IG(#NDIG)
C         C
C               NDA  = #NDA
C               NDU  = #NDU
C               NDF  = #NDF
C               NDIA = #NDIA
C               NDJA = #NDJA
C               NDIG = #NDIG
C         C
C         C     SET UP A, F, IA, JA AND SPECIFY NECESSARY PARAMETERS
C         C
C               ....
C               ....
C         C
C               CALL AMG1R5(A,IA,JA,U,F,IG,
C        +                  NDA,NDIA,NDJA,NDU,NDF,NDIG,NNU,MATRIX,
C        +                  ISWTCH,IOUT,IPRINT,
C        +                LEVELX,IFIRST,NCYC,EPS,MADAPT,NRD,NSOLCO,NRU,
C        +                  ECG1,ECG2,EWT2,NWT,NTR,
C        +                  IERR)
C         C
C               ....
C               ....
C         C
C               STOP
C               END
C
C-----------------------------------------------------------------------
C
C     INPUT VIA ARRAYS (SEE ABOVE):
C
C     A        -   MATRIX L
C
C     IA       -   POINTER VECTOR
C
C     JA       -   POINTER VECTOR
C
C     U        -   FIRST APPROXIMATION TO SOLUTION
C
C     F        -   RIGHT HAND SIDE
C
C
C-----------------------------------------------------------------------
C
C
C     SCALAR INPUT PARAMETERS OF AMG1R5:
C
C     THE INPUT PARAMETERS OF AMG1R5 IN THE LIST BELOW ARE ARRANGED
C     ACCORDING TO THEIR IMPORTANCE TO THE GENERAL USER. THE PARAMETERS
C     PRECEEDED BY A * MUST BE SPECIFIED EXPLICITELY. ALL THE OTHER
C     PARAMETERS ARE SET TO STANDARD VALUES IF ZERO ON INPUT.
C
C     THERE ARE FOUR CLASSES OF INPUT PARAMETERS WITH DECREASING PRI-
C     ORITY:
C
C     1. PARAMETERS DESCRIBING THE USER-DEFINED PROBLEM AND DIMENSIONING
C        OF VECTORS IN THE CALLING PROGRAM
C
C     2. PARAMETERS SPECIFYING SOME GENERAL ALGORITHMIC ALTERNATIVES AND
C        THE AMOUNT OF OUTPUT DURING SOLUTION
C
C     3. PARAMETERS CONTROLLING THE MULTIGRID CYCLING DURING THE SOLU-
C        TION PHASE
C
C     4. PARAMETERS CONTROLLING THE CREATION OF COARSER GRIDS AND INTER-
C        POLATION FORMULAS.
C
C     ONLY THE CLASS 1 - PARAMETERS MUST BE SPECIFIED EXPLICITELY BY
C     THE USER. CLASS 2 - PARAMETERS CONTROL THE GENERAL PERFORMANCE OF
C     AMG1R5. CHANGING THEM DOESN'T REQUIRE UNDERSTANDING THE AMG -
C     ALGORITHM. SPECIFYING NON-STANDARD-VALUES FOR CLASS 3 - PARAMETERS
C     PRESUPPOSES A GENERAL KNOWLEDGE OF MULTIGRID METHODS, WHEREAS THE
C     FUNCTION OF CLASS 4 - PARAMETERS IS ONLY UNDERSTANDABLE AFTER
C     STUDYING THE AMG-ALGORITHM IN DETAIL. FORTUNATELY IN MOST CASES
C     THE CHOICE OF CLASS 3 AND 4 - PARAMETERS ISN'T CRITICAL AND USING
C     THE AMG1R5 - SUPPLIED STANDARD VALUES SHOULD GIVE SATISFACTORY
C     RESULTS.
C
C         --------------------------------------------------------------
C
C     CLASS 1 - PARAMETERS:
C
C  *  NDA      -   DIMENSIONING OF VECTOR A IN CALLING PROGRAM
C
C  *  NDIA     -   DIMENSIONING OF VECTOR IA IN CALLING PROGRAM
C
C  *  NDJA     -   DIMENSIONING OF VECTOR JA IN CALLING PROGRAM
C
C  *  NDU      -   DIMENSIONING OF VECTOR U IN CALLING PROGRAM
C
C  *  NDF      -   DIMENSIONING OF VECTOR F IN CALLING PROGRAM
C
C  *  NDIG     -   DIMENSIONING OF VECTOR IG IN CALLING PROGRAM
C
C  *  NNU      -   NUMBER OF UNKNOWNS
C
C  *  MATRIX   -   INTEGER VALUE CONTAINING INFO ABOUT THE MATRIX L.
C
C                  1ST DIGIT OF MATRIX  --  ISYM:
C                    =1: L IS SYMMETRIC;
C                    =2: L IS NOT SYMMETRIC.
C
C                  2ND DIGIT OF MATRIX  --  IROW0:
C                    =1: L HAS ROWSUM ZERO;
C                    =2: L DOES NOT HAVE ROWSUM ZERO.
C
C         --------------------------------------------------------------
C
C     CLASS 2 - PARAMETERS:
C
C     ISWTCH   -   PARAMETER CONTROLLING WHICH MODULES OF AMG1R5 ARE TO
C                  BE USED.
C                    =1:   CALL FOR -----, -----, -----, WRKCNT.
C                    =2:   CALL FOR -----, -----, SOLVE, WRKCNT.
C                    =3:   CALL FOR -----, FIRST, SOLVE, WRKCNT.
C                    =4:   CALL FOR SETUP, FIRST, SOLVE, WRKCNT.
C                  SETUP DEFINES THE OPERATORS NEEDED IN THE SOLUTION
C                         PHASE.
C                  FIRST INITIALIZES THE SOLUTION VECTOR (SEE PARAMETER
C                         IFIRST).
C                  SOLVE COMPUTES THE SOLUTION BY AMG CYCLING (SEE
C                         PARAMETER NCYC).
C                  WRKCNT PROVIDES THE USER WITH INFORMATION ABOUT
C                         RESIDUALS, STORAGE REQUIREMENTS AND CP-TIMES
C                         (SEE PARAMETER IOUT).
C                  IF AMG1R5 IS CALLED THE FIRST TIME, ISWTCH HAS TO
C                  BE =4. INDEPENDENT OF ISWTCH, SINGLE MODULES CAN BE
C                  BYPASSED BY A PROPER CHOICE OF THE CORRESPONDING
C                  PARAMETER.
C
C     IOUT     -   PARAMETER CONTROLLING THE AMOUNT OF OUTPUT DURING
C                  SOLUTION PHASE:
C
C                  1ST DIGIT: NOT USED; HAS TO BE NON-ZERO.
C
C                  2ND DIGIT:
C                    =0: NO OUTPUT (EXCEPT FOR MESSAGES)
C                    =1: RESIDUAL BEFORE AND AFTER SOLUTION PROCESS
C                    =2: ADD.: STATISTICS ON CP-TIMES AND STORAGE REQUI-
C                        REMENTS
C                    =3: ADD.: RESIDUAL AFTER EACH AMG-CYCLE
C
C     IPRINT   -   PARAMETER SPECIFYING THE FORTRAN UNIT NUMBERS FOR
C                  OUTPUT:
C
C                  1ST DIGIT: NOT USED; HAS TO BE NON-ZERO
C
C                  2ND AND 3RD DIGIT  --  IUP: UNIT NUMBER FOR RESULTS
C
C                  4TH AND 5TH DIGIT  --  IUM: UNIT NUMBER FOR MESSAGES
C
C         --------------------------------------------------------------
C
C     CLASS 3 - PARAMETERS:
C
C     LEVELX   -   MAXIMUM NUMBER OF MG-LEVELS TO BE CREATED (>=1).
C
C     IFIRST   -   PARAMETER FOR FIRST APPROXIMATION.
C
C                  1ST DIGIT OF IFIRST: NOT USED; HAS TO BE NON-ZERO.
C
C                  2ND DIGIT OF IFIRST  --  ITYPU:
C                    =0: NO SETTING OF FIRST APPROXIMATION,
C                    =1: FIRST APPROXIMATION CONSTANT TO ZERO,
C                    =2: FIRST APPROXIMATION CONSTANT TO ONE,
C                    =3: FIRST APPROXIMATION IS RANDOM FUNCTION WITH
C                        THE CONCRETE RANDOM SEQUENCE BEING DETERMINED
C                        BY THE FOLLWING DIGITS.
C
C                  REST OF IFIRST  --  RNDU:
C                    DETERMINES THE CONCRETE RANDOM SEQUENCE USED IN
C                    THE CASE ITYPU=3. (IFIRST=13 IS EQUIVALENT TO
C                    IFIRST=1372815)
C
C     NCYC     -   INTEGER PARAMETER DESCRIBING THE TYPE OF CYCLE TO BE
C                  USED AND THE NUMBER OF CYCLES TO BE PERFORMED.
C
C                  1ST DIGIT OF NCYC  --  IGAM:
C                    =1: V -CYCLE,
C                    =2: V*-CYCLE,
C                    =3: F -CYCLE,
C                    =4: W -CYCLE.
C                  IF NCYC IS NEGATIV, THEN THE APPROXIMATION OF THE
C                  PROBLEM ON THE SECOND FINEST GRID IS COMPUTED BY
C                  IGAM V-CYCLES ON THAT PARTICULAR GRID.
C
C                  2ND DIGIT OF NCYC  --  ICGR:
C                    =0: NO CONJUGATE GRADIENT,
C                    =1: CONJUGATE GRADIENT (ONLY FIRST STEP OF CG),
C                    =2: CONJUGATE GRADIENT (FULL CG).
C
C                  3RD DIGIT OF NCYC  --  ICONV:
C                    CONVERGENCE CRITERION FOR THE USER-DEFINED PROBLEM
C                    (FINEST GRID):
C                    =1: PERFORM A FIXED NUMBER OF CYCLES AS GIVEN BY
C                        NCYCLE (SEE BELOW)
C                    =2: STOP, IF  ||RES|| < EPS
C                    =3: STOP, IF  ||RES|| < EPS * |F|
C                    =4: STOP, IF  ||RES|| < EPS * |U| * |DIAG|
C                    WITH ||RES|| = L2-NORM OF RESIDUAL,
C                           EPS     (SEE INPUT PARAMETER EPS)
C                           |F|   = SUPREMUM NORM OF RIGHT HAND SIDE
C                           |U|   = SUPREMUM NORM OF SOLUTION
C                         |DIAG|  = MAXIMAL DIAGONAL ENTRY IN MATRIX L
C                    NOTE THAT IN ANY CASE THE SOLUTION PROCESS STOPS
C                    AFTER AT MOST NCYCLE CYCLES.
C
C                  REST OF NCYC  --  NCYCLE:
C                    MAXIMAL NUMBER OF CYCLES TO BE PERFORMED (>0) OR
C                    NCYCLE=0: NO CYCLING.
C
C     EPS      -   CONVERGENCE CRITERION FOR SOLUTION PROCESS: (SEE
C                  PARAMETER NCYC). NOTE THAT NO MORE THAN NCYCLE CYCLES
C                  ARE PERFORMED, REGARDLESS OF EPS.
C
C     MADAPT   -   INTEGER VALUE SPECIFYING THE CHOICE OF COARSEST
C                  GRID IN CYCLING:
C
C                  1ST DIGIT OF MADAPT  --  MSEL:
C                    =1: IN CYCLING, ALL GRIDS CONSTRUCTED IN THE SETUP
C                        PHASE ARE USED WITHOUT CHECK.
C                    =2: THE NUMBER OF GRIDS IS AUTOMATICALLY REDUCED
C                        IF THE CONVERGENCE FACTOR ON THE COARSER GRIDS
C                        IS FOUND TO BE LARGER THAN A GIVEN VALUE FAC
C                        (SEE BELOW).
C
C                  REST OF MADAPT  --  FAC
C                        THE REST OF MADAPT DEFINES THE FRACTIONAL PART
C                        OF A REAL NUMBER FAC BETWEEN 0.1 AND 0.99, E.G.
C                        MADAPT=258 MEANS MSEL=2 AND FAC=0.58. IF MADAPT
C                        CONSISTS OF ONLY ONE DIGIT, FAC IS SET TO 0.7
C                        BY DEFAULT.
C
C
C     NRD      -   PARAMETER DESCRIBING RELAXATION (DOWNWARDS):
C
C                  1ST DIGIT OF NRD: NOT USED; HAS TO BE NON-ZERO.
C
C                  2ND DIGIT OF NRD  --  NRDX:
C                    ACTUAL NUMBER OF SMOOTHING STEPS TO BE PERFORMED
C                    THE TYPE OF WHICH IS GIVEN BY THE FOLLOWING DIGITS
C
C                  FOLLOWING DIGITS  --  ARRAY NRDTYP:
C                    =1: RELAXATION OVER THE F-POINTS ONLY
C                    =2: FULL GS SWEEP
C                    =3: RELAXATION OVER THE C-POINTS ONLY
C                    =4: FULL MORE COLOR SWEEP, HIGHEST COLOR FIRST
C
C     NSOLCO   -   PARAMETER CONTROLLING THE SOLUTION ON COARSEST GRID:
C
C                  1ST DIGIT  --  NSC:
C                    =1: GAUSS-SEIDEL METHOD
C                    =2: DIRECT SOLVER (YALE SMP)
C
C                  REST OF NSOLCO  --  NRCX: (ONLY IF NSC=1)
C                  NUMBER OF GS SWEEPS ON COARSEST GRID (>=0).
C                  IF NRCX=0, THEN AS MANY GS SWEEPS ARE PERFORMED
C                  AS ARE NEEDED TO REDUCE THE RESIDUAL BY TWO ORDERS
C                  OF MAGNITUDE. (MAXIMAL 100 RELAXATION SWEEPS)
C
C     NRU      -   PARAMETER FOR RELAXATION (UPWARDS), ANALOGOUS TO NRD.
C
C         --------------------------------------------------------------
C
C     CLASS 4 - PARAMETERS:
C
C     ECG1,ECG2-   REAL PARAMETERS AFFECTING THE CREATION OF COARSER
C     EWT2     -   GRIDS AND/OR THE DEFINITION OF THE INTERPOLATION.
C                  THE CHOICE OF THESE PARAMETERS DEPENDS ON
C                  THE ACTUAL AMG VERSION (SEE SUBROUTINE CRSNG)
C
C     NWT      -   INTEGER PARAMETER AFFECTING THE CREATION OF COARSER
C                  GRIDS AND/OR THE DEFINITION OF THE INTERPOLATION.
C                  THE CHOICE OF THIS PARAMETER DEPENDS ON
C                  THE ACTUAL AMG VERSION (SEE SUBROUTINE CRSNG)
C
C     NTR      -   PARAMETER CONTROLLING COARSE-GRID OPERATOR TRUNCATION
C                    =0: PAIRS OF ZEROES ARE REMOVED FROM COARSE GRID
C                        OPERATORS
C                    =1: NO COARSE-GRID OPERATOR TRUNCATION
C
C
C-----------------------------------------------------------------------
C
C     OUTPUT:
C
C     U        -   CONTAINS THE COMPUTED SOLUTION
C
C
C     IERR     -   ERROR PARAMETER:
C
C                    >0: FATAL ERROR (ABNORMAL TERMINATION OF AMG1R5)
C                    <0: NON-FATAL ERROR (EXECUTION OF AMG1R5 CONTINUES)
C
C                  ERROR CODES IN DETAIL:
C
C                  1. DIMENSIONING TOO SMALL FOR VECTOR
C                        A      (IERR = 1)
C                        IA     (IERR = 2)
C                        JA     (IERR = 3)
C                        U      (IERR = 4)
C                        F      (IERR = 5)
C                        IG     (IERR = 6)
C
C                     NO YALE-SMP BECAUSE OF STORAGE (NDA TOO SMALL):
C                               (IERR = -1)
C                     NO YALE-SMP BECAUSE OF STORAGE (NDJA TOO SMALL):
C                               (IERR = -3)
C                     NO CG BECAUSE OF STORAGE (NDU TOO SMALL):
C                               (IERR = -4)
C                     NO SPACE FOR TRANSPOSE OF INTERPOLATION (NDA OR
C                                                     NDJA TOO SMALL):
C                               (IERR = -1)
C
C                  2. INPUT DATA ERRONEOUS:
C
C                     A-ENTRY MISSING, ISYM = 1:           (IERR = -11)
C                     PARAMETER MATRIX MAY BE ERRONEOUS:   (IERR = -12)
C                     DIAGONAL ELEMENT NOT STORED FIRST:   (IERR =  13)
C                     DIAGONAL ELEMENT NOT POSITIV:        (IERR =  14)
C                     POINTER IA ERRONEOUS:                (IERR =  15)
C                     POINTER JA ERRONEOUS:                (IERR =  16)
C                     PARAMETER ISWTCH ERRONEOUS:          (IERR =  17)
C                     PARAMETER LEVELX ERRONEOUS:          (IERR =  18)
C
C                  3. ERRORS OF THE AMG1R5-SYSTEM (SHOULD NOT OCCUR):
C
C                     TRANSPOSE A-ENTRY MISSING:           (IERR =  21)
C                     INTERPOLATION ENTRY MISSING:         (IERR =  22)
C
C                  4. ALGORITHMIC ERRORS:
C
C                     CG-CORRECTION NOT DEFINED:           (IERR =  31)
C                     NO YALE-SMP BECAUSE OF ERROR IN
C                     FACTORIZATION:                       (IERR = -32)
C
C-----------------------------------------------------------------------
C
C     WORK SPACE:
C
C     THE INTEGER VECTOR IG HAS TO BE PASSED TO AMG1R5 AS WORK SPACE.
C
C-----------------------------------------------------------------------
C
C     DIMENSIONING OF INPUT VECTORS AND WORK SPACE:
C
C     IT'S IMPOSSIBLE TO TELL IN ADVANCE THE EXACT STORAGE REQUIREMENTS
C     OF AMG. THUS, THE FOLLOWING FORMULAS GIVE ONLY REASONABLE GUESSES
C     FOR THE VECTOR LENGTHS WHICH HAVE TO BE DECLARED IN THE CALLING
C     PROGRAM. IN THESE FORMULAS NNA DENOTES THE NUMBER OF NON-ZERO
C     ENTRIES IN THE INPUT-MATRIX L AND NNU IS THE NUMBER OF UNKNOWNS.
C
C     VECTOR         NEEDED LENGTH (GUESS)
C       A               3*NNA + 5*NNU
C       JA              3*NNA + 5*NNU
C       IA              2.2*NNU
C       U               2.2*NNU
C       F               2.2*NNU
C       IG              5.4*NNU
C
C-----------------------------------------------------------------------
C
C
C     STANDARD CHOICES OF PARAMETERS (AS FAR AS MEANINGFUL):
C
C          ISWTCH = 4
C          IOUT   = 12
C          IPRINT = 10606
C
C          LEVELX = 25
C          IFIRST = 13
C          NCYC   = 10110
C          EPS    = 1.D-12
C          MADAPT = 27
C          NRD    = 1131
C          NSOLCO = 110
C          NRU    = 1131
C
C          ECG1   = 0.
C          ECG2   = 0.25
C          EWT2   = 0.35
C          NWT    = 2
C          NTR    = 0
C
C     IF ANY ONE OF THESE PARAMETERS IS 0 ON INPUT, ITS CORRESPONDING
C     STANDARD VALUE IS USED BY AMG1R5.
C
C-----------------------------------------------------------------------
C
C     PORTABILITY RESTRICTIONS:
C
C     1. ROUTINE CTIME IS MACHINE DEPENDENT AND HAS TO BE ADAPTED TO
C        YOUR COMPUTER INSTALLATION OR REPLACED BY A DUMMY ROUTINE.
C
C     2. MOST INPUT PARAMETERS ARE COMPOSED OF SEVERAL DIGITS, THEIR
C        SIGNIFICANCE HAVING BEEN DESCRIBED ABOVE. BE SURE NOT TO ENTER
C        MORE DIGITS THAN YOUR COMPUTER CAN STORE ON AN INTEGER VARI-
C        ABLE.
C
C     3. APART FROM FORTRAN INTRINSIC FUNCTIONS AND SERVICE ROUTINES,
C        THERE IS ONLY ONE EXTERNAL REFERENCE TO A PROGRAM NOT CONTAINED
C        IN THE AMG1R5 - SYSTEM, I.E. THE LINEAR SYSTEM SOLVER NDRV OF
C        THE YALE SPARSE MATRIX PACKAGE. IF YOU HAVN'T ACCESS TO THIS
C        PACKAGE, ENTER A DUMMY ROUTINE NDRV AND AVOID CHOOSING NSC=2
C        (SUBPARAMETER OF NSOLCO). THEN NDRV ISN'T CALLED BY AMG1R5.
C        IN THIS CASE, HOWEVER, INDEFINITE PROBLEMS WILL NOT BE SOLV-
C        ABLE.
C          THE YALE SPARSE MATRIX PACKAGE IS FREELY AVAILABLE FOR NON-
C        PROFIT PURPOSES. CONTACT THE DEPARTMENT OF COMPUTER SCIENCE,
C        YALE UNITVERSITY.
C
C     4. IN AMG1R5 THERE IS THE PARAMETER LRATIO, DENOTING THE RATIO
C        OF SPACE OCCUPIED BY A DOUBLE PRECISION REAL VARIABLE AND
C        THAT OF AN INTEGER. FOR THE IBM-VERSION LRATIO HAS BEEN SET
C        TO 2. CHANGE THIS VALUE IF NECESSARY. (THE SAME HAS TO BE
C        DONE WITH THE YALE SMP-ROUTINE NDRV.)
C
C
C-----------------------------------------------------------------------
C
C     AUTHORS:
C
C          JOHN RUGE, FORT COLLINS (USA),
C              INSTITUTE FOR COMPUTATIONAL STUDIES AT CSU;
C
C          KLAUS STUEBEN, D-5205 ST. AUGUSTIN (W.-GERMANY),
C              GESELLSCHAFT FUER MATHEMATIK UND DATENVERARBEITUNG (GMD).
C
C          ROLF HEMPEL, D-5205 ST. AUGUSTIN (W.-GERMANY),
C              GESELLSCHAFT FUER MATHEMATIK UND DATENVERARBEITUNG (GMD).
C
C-----------------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION A(*),U(*),F(*)
      REAL TIME(20)
      INTEGER IARR(25),IA(*),JA(*),IG(*)
C
C===> LRATIO HAS TO BE SET TO THE NUMBER OF INTEGERS OCCUPYING THE SAME
C     AMOUNT OF STORAGE AS ONE DOUBLE PRECISION REAL.
C
      PARAMETER (LRATIO=2)
C
C===> MAXGR IS THE MAXIMAL NUMBER OF GRIDS. CHANGING THIS UPPER LIMIT
C     JUST REQUIRES CHANGING THE PARAMETER STATEMENT.
C
      PARAMETER (MAXGR=25)
      DOUBLE PRECISION RESI(MAXGR)
      INTEGER IMIN(MAXGR),IMAX(MAXGR),IMINW(MAXGR),IMAXW(MAXGR),
     +        NSTCOL(MAXGR)
      EXTERNAL SETUP,FIRST,SOLVE,WRKCNT
      SAVE
C
      IERR = 0
C
C===> SET PARAMETERS TO STANDARD VALUES, IF NECCESSARY
C
C
      IF (IOUT.NE.0) THEN
        CALL IDEC(IOUT,2,NDIGIT,IARR)
        KOUT = IARR(2)
      ELSE
        KOUT = 2
      ENDIF
C
      IF (ISWTCH.NE.0) THEN
        KSWTCH = ISWTCH
      ELSE
        KSWTCH = 4
      ENDIF
C
      IF (LEVELX.GT.0) THEN
        KEVELX = MIN(LEVELX,MAXGR)
      ELSEIF (LEVELX.LT.0) THEN
        GOTO 70
      ELSE
        KEVELX = MAXGR
      ENDIF
C
      IF (IPRINT.NE.0) THEN
        CALL IDEC(IPRINT,4,NDIGIT,IARR)
        IUP = 10*IARR(2)+IARR(3)
        IUM = IARR(4)
      ELSE
        IUP = 6
        IUM = 6
      ENDIF
      ICGST = NNU+3
      NDICG = (NDIG-ICGST+1)/2
      IF (NDICG.LE.0) GOTO 60
C
      GOTO (10,20,30,40), KSWTCH
      WRITE (IUM,9920)
      IERR = 17
      RETURN
C
40    CALL SETUP (NNU,MATRIX,KEVELX,ECG1,ECG2,EWT2,NWT,NTR,IERR,
     +            A,U,IA,JA,IG,IMIN,IMAX,IMINW,IMAXW,IG(ICGST),
     +            IG(ICGST+NDICG),NSTCOL,IARR,TIME,LEVELS,IROW0,
     +            NDA,NDIA,NDJA,NDU,NDF,NDICG,IUM,MDA,MDIA,
     +            MDJA,MDU,MDF,MDIG,MAXGR,LRATIO)
      IF (IERR.GT.0) RETURN
30    CALL FIRST (IFIRST,U,IMIN,IMAX,IARR,IROW0)
20    CALL SOLVE (MADAPT,NCYC,NRD,NSOLCO,NRU,KOUT,IERR,A,U,F,IA,JA,IG,
     +            EPS,IMIN,IMAX,IMINW,IMAXW,IG(ICGST),IG(ICGST+NDICG),
     +            NSTCOL,IARR,TIME,NCYC0,IROW0,LEVELS,NDA,NDJA,
     +            NDU,NDF,MDA,MDJA,MDU,MDF,IUP,IUM,RESI,RES0,RES)
      IF (IERR.GT.0) RETURN
10    CALL WRKCNT(KOUT,IA,IG,IMIN,IMAX,IMINW,LEVELS,TIME,NCYC0,IUP,
     +            MDA,MDIA,MDJA,MDU,MDF,MDIG,RES0,RES)
      RETURN
60    WRITE (IUM,9910)
      IERR = 6
      RETURN
70    WRITE (IUM,9930)
      IERR = 18
      RETURN
9910  FORMAT (' *** ERROR IN AMG1R5: NDIG TOO SMALL ***')
9920  FORMAT (' *** ERROR IN AMG1R5: ILLEGAL PARAMETER ISWTCH ***')
9930  FORMAT (' *** ERROR IN AMG1R5: ILLEGAL PARAMETER LEVELX ***')
      END
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     GENERAL AMG1R5 SETUP-SUBROUTINES
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C.......................................................................
C
C     SETUP                                                  SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE SETUP(NNU,MATRIX,LEVELX,ECG1,ECG2,EWT2,NWT,NTR,IERR,
     +                 A,U,IA,JA,IW,IMIN,IMAX,IMINW,IMAXW,ICG,IFG,
     +                 NSTCOL,IARR,TIME,LEVELS,IROW0,NDA,NDIA,NDJA,
     +                 NDU,NDF,NDICG,IUM,MDA,MDIA,MDJA,MDU,MDF,MDIG,
     +                 MAXGR,LRATIO)
C
C     PREPARATION PHASE OF AMG1R5 (GENERAL PART)
C
      DOUBLE PRECISION A(*),U(*),ECG1,ECG2,EWT2
      REAL TIME(*)
      INTEGER IARR(*),IA(*),JA(*),IW(*),IMIN(*),IMAX(*),IMINW(*),
     +        IMAXW(*),ICG(*),IFG(*),NSTCOL(*)
      EXTERNAL IDEC,CHECK,CRSNG
      SAVE
C
C===> DECOMPOSE "MATRIX"
C
      CALL IDEC(MATRIX,2,NDIGIT,IARR)
      ISYM  = IARR(1)
      IROW0 = IARR(2)
C
C===> PREPARATION (IGNORED IN TIMING)
C
      IMIN(1) = 1
      IMAX(1) = NNU
      CALL CHECK(IERR,A,IA,JA,IMIN,IMAX,ICG,IFG,TIME,ISYM,IROW0,NDA,
     +           NDICG,NDJA,IUM)
      IF (IERR.GT.0) RETURN
C
C===> RESET TIME COUNTERS
C
      DO 30 I=1,20
        TIME(I) = 0.0
30    CONTINUE
C
C===> DEFINE COARSER GRIDS + OPERATORS. RESET LEVELS IF NECESSARY.
C
      CALL CRSNG(LEVELX,ECG1,ECG2,EWT2,NWT,NTR,IERR,A,U,IA,JA,IW,IMIN,
     +           IMAX,IMINW,IMAXW,ICG,IFG,NSTCOL,LEVELS,IROW0,NDA,
     +           NDJA,NDIA,NDU,NDF,NDICG,TIME,IUM,MDA,MDIA,MDJA,MDU,MDF,
     +           MDIG,MAXGR,LRATIO)
      RETURN
      END
C
C.......................................................................
C
C     CHECK                                                 SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE CHECK(IERR,A,IA,JA,IMIN,IMAX,ICG,IFG,TIME,ISYM,IROW0,
     +                 NDA,NDICG,NDJA,IUM)
C
C     CHECKS FOR SEVERAL PROPERTIES OF A, IA, JA. IN PARTICULAR,
C     CHECKS FOR SYMMETRIC STORAGE OF GIVEN MATRIX (I.E. L(I,J) IS
C     STORED IN A IFF L(J,I) IS STORED IN A).
C
C     IF STORAGE IS SYMMETRIC, PAIRS OF ZEROES ARE REMOVED (IF THERE
C     ARE ANY) AND PROGRAM EXECUTION CONTINUES.
C
C     IF, HOWEVER, STORAGE OF A IS NOT SYMMETRIC, IT IS SYMMETRIZED.
C     IF ISYM.EQ.1, THE MISSING TRANSPOSE CONNECTIONS ARE COPIED.
C     IF ISYM.NE.1, THE STORAGE OF THE MATRIX A IS SYMMETRIZED (BY
C     ADDING CERTAIN ZERO ELEMENTS). THEN PAIRS OF ZEROES ARE REMOVED
C     (IF THERE ARE ANY) AND PROGRAM EXECUTION CONTINUES.
C
C     ARRAYS USED FOR TEMPORARY STORAGE: ICG, IFG
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION A(*)
      REAL TIME(*)
      INTEGER IA(*),JA(*),IMIN(*),IMAX(*),ICG(*),IFG(*)
      EXTERNAL TRUNC
      SAVE
C
C===> CHECK IF POINTERS IA, JA ARE REASONABLE
C
      IF (IA(1).NE.1) GOTO 3300
      NNU = IMAX(1)
      IF (NNU.GE.NDICG) GOTO 3200
      DO 2 I=1,NNU
        ICG(I) = 0
2     CONTINUE
      DO 10 I=1,NNU
        J1 = IA(I)
        J2 = IA(I+1)-1
        IF (J2.LT.J1 .OR. J2.GT.NDA) GOTO 3300
        IF (JA(J1).NE.I)             GOTO 3320
        DO 5 J=J1,J2
          I1 = JA(J)
          IF (I1.LT.1 .OR. I1.GT.NNU .OR. ICG(I1).EQ.1) GOTO 3310
          ICG(I1) = 1
5       CONTINUE
        DO 7 J=J1,J2
          ICG(JA(J)) = 0
7       CONTINUE
10    CONTINUE
      NNA = IA(NNU+1)-1
C
C===> CHECK FOR PROPERTIES OF A. IN PARTICULAR, COUNT
C===> MISSING STORAGE PLACES ("NEW"). RETURN IF NEW=0
C
      ANORMM = 0.0D0
      ANORMP = 0.0D0
      NAOFF  = 0
      NAPOS  = 0
      NANEG  = 0
      NAZER  = 0
      NEW    = 0
      DO 100 I=1,NNU
         ICG(I) = IA(I+1)-IA(I)
100   CONTINUE
C
      DO 200 I=1,NNU
        D = A(IA(I))
        IF (D.LE.0.0D0) GOTO 3330
        DEPS = D*1.0D-12
        ROWSUM = D
        ANORMP = ANORMP+2.0D0*D**2
        DO 150 J=IA(I)+1,IA(I+1)-1
          ROWSUM = ROWSUM+A(J)
          IF (A(J).GE.DEPS) NAOFF = NAOFF+1
          I1 = JA(J)
          IF (I1.GT.0) GOTO 140
          JA(J) = -I1
          GOTO 150
140       IF(I1.LT.I) GOTO 135
          DO 130 J1=IA(I1)+1,IA(I1+1)-1
            IF (JA(J1).NE.I) GOTO 130
            JA(J1) = -JA(J1)
            ANORMM = ANORMM+(A(J)-A(J1))**2
            ANORMP = ANORMP+(A(J)+A(J1))**2
            GOTO 150
130       CONTINUE
135       JA(J)=-JA(J)
          ANORMM = ANORMM+A(J)**2
          ANORMP = ANORMP+A(J)**2
          NEW = NEW+1
          ICG(I1) = ICG(I1)+1
150     CONTINUE
        IF (ROWSUM.GT.DEPS) THEN
          NAPOS = NAPOS+1
        ELSE IF (ROWSUM.LT.-DEPS) THEN
          NANEG = NANEG+1
        ELSE
          NAZER = NAZER+1
        ENDIF
200   CONTINUE
      ANORMM = SQRT(ANORMM)
      ANORMP = SQRT(ANORMP)
      ASYM = ANORMM/ANORMP
      IF (ASYM.LE.1.0D-12) ASYM = 0.0D0
C
C===> MESSAGES ON A
C
      IF (ASYM.EQ.0.0D0) WRITE (IUM,9000)
      IF (ASYM.NE.0.0D0) WRITE (IUM,9005) ASYM
      IF (NAOFF.GT.0)    WRITE (IUM,9010) NAOFF
      IF (NANEG.GT.0)    WRITE (IUM,9020) NANEG
      IF (NAZER.EQ.NNU)  WRITE (IUM,9025)
      IF (NAOFF.EQ.0 .AND. NANEG.EQ.0) WRITE (IUM,9030)
C
C===> WARNINGS
C
      IF (ISYM.EQ.1  .AND. ASYM.NE.0.0D0 .OR.
     +    ISYM.GT.1  .AND. ASYM.EQ.0.0D0 .OR.
     +    IROW0.EQ.1 .AND. NAZER.NE.NNU  .OR.
     +    IROW0.GT.1 .AND. NAZER.EQ.NNU) THEN
        WRITE (IUM,9100)
        IERR = -12
      ENDIF
C
      IF (NEW.GT.0) GOTO 220
      WRITE (IUM,9600)
      GOTO 600
C
C===> REPLACE A BY SYMMETRIZED VERSION
C
220   IF (NNA+NEW.GE.NDA) GOTO 3210
      IF (NNA+NEW.GE.NDJA) GOTO 3220
C
C===> EXTEND MATRIX A IN SITU
C
      IFG(1) = 1
      DO 230 I=2,NNU+1
        IFG(I) = IFG(I-1)+ICG(I-1)
230   CONTINUE
      DO 250 I=NNU,1,-1
        ISHIFT = IFG(I)-IA(I)
        DO 240 J=IA(I+1)-1,IA(I),-1
          A(J+ISHIFT)  =  A(J)
          JA(J+ISHIFT) = JA(J)
240     CONTINUE
        ICG(I)  = IA(I+1)+ISHIFT
        IA(I+1) = IFG(I+1)
250   CONTINUE
C
C===> SYMMETRIZE MATRIX A:  ISYM=1: COPY MISSING TRANSPOSE ENTRIES
C                           ISYM=2: FILL IN ZEROES
C
      IF (ISYM.NE.1) THEN
        DO 400 I=1,NNU
          J2 = ICG(I)-1
          DO 450 J=IA(I)+1,J2
            I1 = JA(J)
            IF (I1.GT.0) GOTO 450
            JA(J) = -I1
            JNEW=ICG(JA(J))
            A(JNEW)= 0.0D0
            JA(JNEW) = I
            ICG(JA(J))=JNEW+1
450       CONTINUE
400     CONTINUE
        WRITE (IUM,9500) NEW
      ELSE
        WRITE (IUM,9220) NEW
        IERR = -11
        DO 500 I=1,NNU
          J2 = ICG(I)-1
          DO 550 J=IA(I)+1,J2
            I1 = JA(J)
            IF (I1.GT.0) GOTO 550
            JA(J) = -I1
            JNEW=ICG(JA(J))
            A(JNEW)= 0.0D0
            JA(JNEW) = I
            ICG(JA(J))=JNEW+1
550       CONTINUE
500     CONTINUE
      ENDIF
C
C===> REMOVE PAIRS OF ZEROES
C
600   CALL TRUNC(1,0,A,IA,JA,IMIN,IMAX,TIME,IERR,IUM)
      RETURN
C
C===> ERROR MESSAGES
C
3200  WRITE (IUM,9200)
      IERR = 6
      RETURN
C
3210  WRITE (IUM,9210)
      IERR = 1
      RETURN
C
3220  WRITE (IUM,9230)
      IERR = 3
      RETURN
C
3300  WRITE (IUM,9300)
      IERR = 15
      RETURN
C
3310  WRITE (IUM,9310)
      IERR = 16
      RETURN
C
3320  WRITE (IUM,9320)
      IERR = 13
      RETURN
C
3330  WRITE (IUM,9330)
      IERR = 14
      RETURN
C
9000  FORMAT (' CHECK: A PROBABLY SYMMETRIC')
9005  FORMAT (' CHECK: A PROBABLY NOT SYMMETRIC. MEASURE:',D11.3)
9010  FORMAT (' CHECK: A PROBABLY NOT POS. TYPE:',I6,' OFF-DIAGONAL',
     +        ' ELEMENTS POSITIVE')
9020  FORMAT (' CHECK: A PROBABLY NOT POS. TYPE:',I6,' ROWSUMS',
     +        ' NEGATIVE')
9025  FORMAT (' CHECK: A PROBABLY SINGULAR - ROWSUMS ARE ZERO')
9030  FORMAT (' CHECK: A PROBABLY POSITIVE TYPE')
9100  FORMAT (' --- WARNG IN CHECK: PARAM MATRIX MAY BE BAD ---')
C
9200  FORMAT (' *** ERROR IN CHECK: NDIG TOO SMALL ***')
9210  FORMAT (' *** ERROR IN CHECK: NDA TOO SMALL ***')
9220  FORMAT (' *** WARNG IN CHECK:',I5,' A-ENTRIES MISSING ***'/
     +        '     MISSING TRANSPOSE CONNECTIONS WILL BE FILLED IN')
9230  FORMAT (' *** ERROR IN CHECK: NDJA TOO SMALL ***')
9300  FORMAT (' *** ERROR IN CHECK: POINTER IA ERRONEOUS ***')
9310  FORMAT (' *** ERROR IN CHECK: POINTER JA ERRONEOUS ***')
9320  FORMAT (' *** ERROR IN CHECK: DIAGONAL IS NOT STORED FIRST ***')
9330  FORMAT (' *** ERROR IN CHECK: DIAGONAL IS NON-POSITIVE ***')
9500  FORMAT (' CHECK: STORAGE OF A HAS BEEN SYMMETRIZED BY',
     +        ' INTRODUCING',I6,' ZEROES')
9600  FORMAT (' CHECK: MATRIX A WAS SYMMETRICALLY STORED')
      END
C
C.......................................................................
C
C     TRUNC                                                 SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE TRUNC(K,NTR,A,IA,JA,IMIN,IMAX,TIME,IERR,IUM)
C
C     TRUNCATES OPERATOR ON GRID K CORRESPONDING TO THE VALUE OF NTR.
C     NTR HAS TO BE 0 OR 1:
C
C       =0:    PAIRS OF ZEROES ARE REMOVED FROM COARSE GRID OPERATORS;
C       =1:    COARSE GRID OPERATORS REMAIN UNCHANGED.
C
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION A(*)
      REAL TOLD,TNEW,TIME(*)
      INTEGER IA(*),JA(*),IMIN(*),IMAX(*)
      EXTERNAL CTIME
      SAVE
C
      IF (NTR.EQ.1) RETURN
C
      CALL CTIME(TOLD)
      IMN = IMIN(K)
      IMX = IMAX(K)
      NNA = IA(IMX+1)-IA(IMN)
      JPOS = IA(IMN)
C
      DO 205 I=IMN,IMX
        J1 = IA(I)+1
        J2 = IA(I+1)-1
        A(JPOS)  = A(IA(I))
        JA(JPOS) = I
        IA(I) = JPOS
        JPOS = JPOS+1
        DO 250 J=J1,J2
          I1 = JA(J)
          IF (I1.LT.0) GOTO 250
          IF (I1.LT.I) GOTO 230
          IF (A(J).NE.0.0D0) GOTO 230
          DO 210 J3=IA(I1),IA(I1+1)-1
            IF (JA(J3).NE.I) GOTO 210
            JT = J3
            AT = A(J3)
            GOTO 215
210       CONTINUE
          GOTO 1000
215       IF (AT.NE.0.0D0) GOTO 230
          JA(JT) = -JA(JT)
          GOTO 250
230       A(JPOS)  = A(J)
          JA(JPOS) = JA(J)
          JPOS = JPOS+1
250     CONTINUE
205   CONTINUE
      IA(IMX+1) = JPOS
C
C===> EXIT
C
      CALL CTIME(TNEW)
      TIME(7) = TIME(7)+TNEW-TOLD
      RETURN
C
C===> ERROR MESSAGE
C
1000  WRITE (IUM,9000) K
      IERR = 21
      RETURN
C
9000  FORMAT (' *** ERROR IN TRUNC: TRANSPOSE A-ENTRY MISSING ON GRID'
     +        ,I3,' ***')
      END
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     AMG1R5 SETUP ROUTINES (VERSION 3)
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C.......................................................................
C
C     CRSNG                                           SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE CRSNG(LEVELX,EECG1,EECG2,EEWT2,NNWT,NTR,IERR,A,U,
     +                 IA,JA,IW,IMIN,IMAX,IMINW,IMAXW,ICG,IFG,NSTCOL,
     +                 LEVELS,IROW0,NDA,NDJA,NDIA,NDU,NDF,NDICG,TIME,
     +                 IUM,MDA,MDIA,MDJA,MDU,MDF,MDIG,MAXGR,LRATIO)
C
C     PERFORMS COARSENING, DEFINES INTERPOLATIONS AND CGC-OPERATORS
C
C     ============== STANDARD VALUES OF PARAMETERS =====================
C
C             ECG1=0.,    ECG2=0.25,    EWT2=0.35,   NWT=2
C
C     ================ DESCRIPTION OF PARAMETERS =======================
C
C     ECG1 --    DEFINES CRITERION FOR DETERMINING DIAGONAL
C                DOMINANCE IN RWSRT. I.E., IF THE ABSOLUTE
C                VALUE OF THE SUM OF THE OFF-DIAGONALS OF ROW
C                I IS SMALLER THAN ECG1 TIMES THE ABSOLUTE
C                VALUE OF THE DIAGONAL ENTRY, THEN POINT I IS
C                IMMEDIATELY FORCED TO BE AN F-POINT IN THE
C                PRE-COLORING ALGORITHM PCOL.
C                IN THE SECOND PART (WINT), NO INTERPOLATION
C                IS DEFINED FOR POINT I, AND NO POINTS USE
C                I FOR INTERPOLATION. IN ADDITION, THE WEIGHT
C                FOR POINT I IS NOT DISTRIBUTED TO OTHER POINTS
C                WHEN DEFINING THE INTERPOLATION WEIGHTS FOR
C                POINTS WHICH DEPEND ON POINT I. (THIS IS
C                EQUIVALENT TO COMPLETELY IGNORING SUCH POINTS
C                IN DETERMINING THE COARSE GRID AND INTERPOLATION
C                WEIGHTS.
C
C     ECG2 --    DEFINES STRONG CONNECTIONS (ALPHA IN THE PAPER)
C
C     EWT2 --    DEFINES STRONG DEPENDENCE ON A SET (BETA IN THE PAPER)
C
C     NWT  --    PARAMETER CONTROLLING THE DEFINITION OF INTERPOLATION
C                FORMULAS:
C                  =1 - CHECKING OF INT-FORMULA: OFF
C                  =2 - CHECKING OF INT-FORMULA: ON
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION A(*),U(*)
      REAL TIME(*)
      INTEGER IA(*),JA(*),IW(*),IMIN(*),IMAX(*),IMINW(*),IMAXW(*),
     +        ICG(*),IFG(*),NSTCOL(*)
      EXTERNAL RWSRT,PCOL,PWINT,OPDFN,TRUNC,SETIFG
      SAVE
C
C===> ASSIGN DEFAULT VALUES IF ZERO
C
      NWT = NNWT
      ECG1 = EECG1
      ECG2 = EECG2
      EWT2 = EEWT2
      IF (NWT.EQ.0)     NWT = 2
      IF (ECG2.EQ.0.0D0) ECG2 = 0.25D0
      IF (EWT2.EQ.0.0D0) EWT2 = 0.35D0
C
C===> DECODE PARAMETER NWT
C
    5 ICHK = NWT
C
C===> COARSENING
C
   10 LEVELS = MIN(LEVELX,MAXGR)
      MMAX   = MAXGR
C
C===> INITIALIZE PARAMETERS MDA - MDIW, LATER SET TO ACTUAL STORAGE
C     REQUIREMENTS OF CORRESPONDING VECTORS
C
      MDA   = 0
      MDJA  = 0
      MDICG = 0
      MDIR  = 0
      MDIW  = IMAX(1)
C
C===> KFIRST IS THE NUMBER OF THE FIRST STORED GRID, IAJAS, IIAS AND
C     IIRS ARE THE SHIFTS IN VECTORS A, JA, IA AND IR, RESPECTIVELY, AS
C     COMPARED TO FULL STORAGE OF ALL GRIDS
C
      KERR   = 0
      KFIRST = 1
      IAJAS  = 0
      IIAS   = 0
      IIRS   = 0
C
      DO 50 K=2,LEVELS
C
C===> JTRST: INITIAL POINTER FOR WORK SPACE IN VECTOR A, TO CONTAIN
C            THE STRONG TRANSPOSE CONNECTIONS
C     NDJTR: AVAILABLE WORK SPACE
C
20      JTRST = IA(IMAX(K-1)+1)
        NDJTR = LRATIO*(NDA-JTRST+1)
        CALL RWSRT(K-1,ECG1,ECG2,IERR,A,IA,JA,IW,IMIN,IMAX,IMAXW,
     +             IFG,A(JTRST),TIME,NDJTR,IUM,MDJTR)
        IF (IERR.GT.0) THEN
          IF (IERR.GE.1.AND.IERR.LE.6) GOTO 30
          RETURN
        ENDIF
        MDA = MAX(MDA,JTRST+(MDJTR-1)/LRATIO+IAJAS)
C
C===>   IRST: INITIAL POINTER FOR WORK SPACE IN VECTOR U, TO CONTAIN
C             RESET STACK
C       NDIR: AVAILABLE WORK SPACE
C
        IRST = MDIW+1-IIRS
        NDIR = LRATIO*(NDU-IRST+1)
        CALL PCOL(K-1,IERR,IA,JA,IW,IMIN,IMAX,IMAXW,ICG,IFG,U(IRST),
     +            A(JTRST),TIME,NDICG,NDIR,IUM,MDICG,MDIR,IIAS)
        IF (IERR.GT.0) THEN
          IF (IERR.GE.1.AND.IERR.LE.6) GOTO 30
          RETURN
        ENDIF
        CALL PWINT(K-1,EWT2,ICHK,MMAX,A,IA,JA,IW,IMIN,IMAX,IMINW,
     +             IMAXW,IFG,ICG,U(IRST),TIME,IERR,IROW0,NCOLX,
     +             NDA,NDJA,NDICG,NDIA,IUM,MDA,MDJA,IAJAS)
        IF (IERR.GT.0) THEN
          IF (IERR.GE.1.AND.IERR.LE.6) GOTO 30
          RETURN
        ENDIF
C
        CALL OPDFN(K,IERR,MMAX,A,IA,JA,IW,IMIN,IMAX,IMINW,IMAXW,ICG,
     +             IFG,U(IRST),NSTCOL,NCOLX,TIME,NDA,NDJA,IUM,
     +             MDA,MDJA,IAJAS)
        IF (IERR.GT.0) THEN
          IF (IERR.GE.1.AND.IERR.LE.6) GOTO 30
          RETURN
        ENDIF
        IF (K.GT.MMAX) GOTO 100
        CALL TRUNC(K,NTR,A,IA,JA,IMIN,IMAX,TIME,IERR,IUM)
        IF (IERR.GT.0) RETURN
        IW(IMINW(K)) = IA(IMAX(K)+1)
        IF (K.GE.MMAX) GOTO 100
        GOTO 50
30      IF (K.LE.KFIRST+1.AND.IIRS.NE.0) RETURN
        KERR = IERR
        IERR = 0
        KFIRST = K-1
        IIRS = MDIW
        ISIA = IMIN(K-1)-1
        IIAS = IIAS+ISIA
        ISAJA = IA(IMIN(K-1))-1
        IAJAS = IAJAS+ISAJA
        DO 35 I=IA(IMIN(K-1)),IA(IMAX(K-1)+1)-1
          JA(I-ISAJA) = JA(I)-ISIA
          A(I-ISAJA)  = A(I)
35      CONTINUE
        DO 40 I=IMIN(K-1),IMAX(K-1)+1
          IA(I-ISIA) = IA(I)-ISAJA
40      CONTINUE
        IMIN(K-1) = 1
        IMAX(K-1) = IMAX(K-1)-ISIA
        IW(IMINW(K-1)) = IA(IMAX(K-1)+1)
        GOTO 20
50    CONTINUE
      GOTO 200
100   LEVELS=MMAX
200   MDIA = IMAX(LEVELS)+1+IIAS
      MDU  = MAX(IMAX(LEVELS)+IIAS,MDIW+1+(MDIR-1)/LRATIO,
     +                         MDIW+1+(MDIW-1)/LRATIO)
      MDF  = IMAX(LEVELS)+IIAS
      MDIG = MDIW+2+2*MDICG
      IF (KERR.NE.0.OR.MDU.GT.NDU.OR.MDF.GT.NDF) THEN
        WRITE (IUM,9000) MDA,MDJA,MDIA,MDU,MDF,MDIG
        WRITE (IUM,9010) MDIW
        IERR = KERR
        RETURN
      ENDIF
      CALL SETIFG(IMIN,IMAX,ICG,IFG,NSTCOL,LEVELS,TIME)
      RETURN
9000  FORMAT (/' **************** SPACE REQUIREMENTS ****************'//
     +         ' VECTOR          NEEDED                              '/
     +         ' ----------------------------------------------------'/
     +         '    A ',       I16,   '   ADJUST THE DIMENSIONING OF '/
     +         '    JA',       I16,   '   VECTORS A - IG IN THE      '/
     +         '    IA',       I16,   '   CALLING PROGRAM ACCORDING  '/
     +         '    U ',       I16,   '   TO THE CALCULATED SPACE RE-'/
     +         '    F ',       I16,   '   QUIREMENTS AND RERUN THE   '/
     +         '    IG',       I16,   '   PROGRAM.                   '/
     +         ' ----------------------------------------------------')
9010  FORMAT (/' NOTE: IF YOU WANT TO USE CG-CORRECTIONS IN THE SOLU-'/
     +         '       TION PROCESS (NCYC-SUBPARAMETER ICGR=1 OR =2),'/
     +         '       PROVIDE FOR ADDITIONAL',I6,' STORAGE LOCATIONS'/
     +         '       IN VECTORS U AND F.                           '/
     +         '         SIMILARLY, USAGE OF THE YALE-SMP SOLVER ON  '/
     +         '       THE COARSEST GRID (NSOLCO=2) WILL REQUIRE     '/
     +         '       ADDITIONAL SPACE IN VECTOR A DURING THE SOLU- '/
     +         '       TION PHASE. IN THIS CASE, HOWEVER, ITS EXACT  '/
     +         '       AMOUNT ISN''T PREDICTABLE.                    '/)
      END
C
C.......................................................................
C
C     RWSRT                                             SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE RWSRT(K,ECG1,ECG2,IERR,A,IA,JA,IW,IMIN,IMAX,
     +                 IMAXW,IFG,JTR,TIME,NDJTR,IUM,MDJTR)
C
C     ROW-SORT ALGORITHM FOR ROWS OF A(K). IN DETAIL:
C
C     - SORTS THE ELEMENTS OF EACH ROW OF A(K) SUCH THAT THE STRONG
C       CONNECTIONS JUST FOLLOW THE DIAGONAL ELEMENT. ONE OF THE STRONG-
C       EST IS ALWAYS FIRST (EVEN IF PARTIAL SORTING IS PERFORMED!).
C       NON-NEGATIVE CONNECTIONS ARE ALWAYS DEFINED TO BE WEAK.
C       JA(IA(I)) IS RE-DEFINED TO POINT TO THE LAST STRONG CONNECTION
C       OF POINT I (OR TO IA(I) IF THERE IS NO STRONG CONNECTION).
C
C     - THE STRONG TRANSPOSE CONNECTIONS ARE LOADED LOGICALLY
C       INTO JTR, I.E., ITS POINTERS ARE STORED IN JTR.
C
C     ============== COMMENTS ON WORK SPACE USED =======================
C
C     - JA(IA(I)) ----- I=IMIN(K),...,IMAX(K)
C
C     IS DEFINED TO POINT TO THE LAST STRONG CONNECTION OF POINT I.
C     NOTE THAT THE ORIGINAL CONTENTS OF JA(IA(I)) IS OVERWRITTEN.
C     (IT IS PUT BACK IN SUBROUTINE WINT.)
C
C     - JTR(J)--------- J=1,IW(IMAX(K)+IWS+1)-1
C     - IW(I) --------- I=IMIN(K)+IWS,...,IMAX(K)+IWS
C     - IFG(I)--------- I=IMIN(K),...,IMAX(K)
C
C     (IWS=0, IF K=1; IWS=IMAXW(K-1)+2-IMIN(K) OTHERWISE)
C
C     JTR IS INITIALIZED TO HAVE SAME FORM AS JA. JTR(J)
C     CONTAINS INFORMATION ON STRONG TRANSPOSE CONNECTIONS:
C     JTR(J) WITH IW(I+IWS)<=J<=IW(I+IWS+1)-1 POINTS TO THE STRONG
C     TRANSPOSE CONNECTIONS OF I.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION A(*)
      REAL TOLD,TNEW,TIME(*)
      INTEGER IA(*),JA(*),IW(*),JTR(*),IMIN(*),IMAX(*),IMAXW(*),IFG(*)
      EXTERNAL CTIME
      SAVE
C
C===> INITIALIZATION OF WORK SPACE
C
      CALL CTIME(TOLD)
      ILO = IMIN(K)
      IHI = IMAX(K)
      IF (K.NE.1) THEN
        IWS = IMAXW(K-1)+2-ILO
      ELSE
        IWS = 0
      ENDIF
      DO 5 I=ILO,IHI
        IFG(I) = 0
        JA(IA(I)) = IA(I)
    5 CONTINUE
      CALL CTIME(TNEW)
      TIME(8) = TIME(8)+TNEW-TOLD
      TOLD = TNEW
C
C     *************************
C     * PARTIAL STANDARD SORT *
C     *************************
C
C     NOTE: NON-NEGATIVE CONNECTIONS ARE ALWAYS DEFINED TO BE WEAK! THE
C           STRONGEST CONNECTION IS ALWAYS FOLLOWING THE DIAGONAL.
C
      DO 590 I=ILO,IHI
        JLO = IA(I)+1
        JHI = IA(I+1)-1
        IF (JHI.LT.JLO) GOTO 590
C
C===>   FIND STRONGEST CONNECTION AND SUM OF |OFF-DIAGONALS| OF ROW I
C
        AMX = A(JLO)
        AMN = A(JLO)
        JMX = JLO
        IF (ECG1.NE.0.D0) THEN
          RS = 0.0D0
          DO 550 J=JLO+1,JHI
            RS = RS+ABS(A(J))
            IF (A(J).LT.AMX) THEN
              AMX = A(J)
              JMX = J
            ELSEIF (A(J).GT.AMN) THEN
              AMN = A(J)
            ENDIF
  550     CONTINUE
C
C===>     TEST FOR POSITIVE OFF-DIAGONALS / DIAGONAL DOMINANCE
C
          IF (AMX.GE.0.0D0 .OR. RS.LE.ECG1*A(IA(I))) GOTO 590
        ELSE
          DO 555 J=JLO+1,JHI
            IF (A(J).LT.AMX) THEN
              AMX = A(J)
              JMX = J
            ELSEIF (A(J).GT.AMN) THEN
              AMN = A(J)
            ENDIF
  555     CONTINUE
C
C===>   TEST FOR POSITIVE OFF-DIAGONALS
C
          IF (AMX.GE.0.D0) GOTO 590
        ENDIF
C
C===>   PUT STRONGEST CONNECTION IN FIRST POSITION
C
        AST = ECG2*AMX
        IMX = JA(JMX)
        A(JMX)  = A(JLO)
        JA(JMX) = JA(JLO)
        A(JLO)  = AMX
        JA(JLO) = IMX
        IF (AMN.LE.AST) GOTO 580
        JHI = JHI+1
C
C===>   DECREASE JHI UNTIL A STRONG CONNECTION IS FOUND
C       (IF JLO >= JHI STOP: ALL CONNECTIONS ARE SORTED)
C
  560   JHI = JHI-1
        IF (JLO.GE.JHI) GOTO 580
        IF (A(JHI).GT.AST) GOTO 560
C
C===>   INCREASE JLO UNTIL A WEAK CONNECTION IS FOUND
C       (IF JLO >= JHI STOP: ALL CONNECTIONS ARE SORTED)
C
  570   JLO = JLO+1
        IF (JLO.GE.JHI) GOTO 580
        IF (A(JLO).LE.AST) GOTO 570
C
C===>   INTERCHANGE A(JHI) AND A(JLO)
C
        ATMP = A(JHI)
        ITMP = JA(JHI)
        A(JHI)  = A(JLO)
        JA(JHI) = JA(JLO)
        A(JLO)  = ATMP
        JA(JLO) = ITMP
        GOTO 560
C
C===>   ROW SORTED --  SET JA(IA(I)) TO LAST STRONG CONNECTION
C
  580   JA(IA(I)) = JHI
C
C===>   COUNT STRONG TRANSPOSE CONNECTIONS IN ROW I
C
        DO 585 J=IA(I)+1,JHI
          IFG(JA(J)) = IFG(JA(J))+1
  585   CONTINUE
  590 CONTINUE
C
C===> INITIALIZATION OF WORK SPACE FOR STRONG TRANSPOSE CONNECTIONS
C
      IW(ILO+IWS) = 1
      DO 5010 I=ILO,IHI
        IW(I+IWS+1) = IW(I+IWS)+IFG(I)
        IFG(I) = IW(I+IWS)
 5010 CONTINUE
      MDJTR = IW(IHI+IWS+1)-1
      IF (MDJTR.GT.NDJTR) GOTO 9901
C
C===> LOAD POINTERS TO STRONG TRANSPOSE CONNECTIONS INTO JTR
C
      DO 5030 I=ILO,IHI
        DO 5020 J=IA(I)+1,JA(IA(I))
          II = JA(J)
          JTR(IFG(II)) = I
          IFG(II) = IFG(II)+1
 5020   CONTINUE
 5030 CONTINUE
      CALL CTIME(TNEW)
      TIME(1) = TIME(1)+TNEW-TOLD
      RETURN
C
C===> ERROR MESSAGES
C
 9901 WRITE (IUM,9910)
      IERR = 1
      RETURN
C
 9910 FORMAT (' *** ERROR IN RWSRT: NDA TOO SMALL ***')
      END
C
C.......................................................................
C
C     PCOL                                            SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE PCOL(K,IERR,IA,JA,IW,IMIN,IMAX,IMAXW,ICG,IFG,IR,JTR,
     +                TIME,NDICG,NDIR,IUM,MDICG,MDIR,IIAS)
C
C     PRE-COLORING ALGORITHM FOR GRID K. THIS IS THE VERSION AS
C     DESCRIBED IN RUGE/STUEBEN (BRISTOL). THE GOAL IS TO OBTAIN QUICKLY
C     A TENTATIVE SET OF C-POINTS WITH THE FOLLOWING PROPERTIES:
C
C       - THE C-POINTS ARE ONLY WEAKLY CONNECTED AMONG EACH OTHER;
C       - EACH F-POINT HAS (AT LEAST) ONE STRONG CONNECTION TO A C-POINT
C         (EXCEPT FOR SOME EXCEPTIONAL F-POINTS, E.G., THOSE WHICH DO
C         NOT HAVE A STRONG CONNECTION AT ALL: THE "FORCED F-POINTS").
C
C     ON EXIT, ICG(I) (IMIN(K)<=I<=IMAX(K)) CONTAINS ALL THE INFORMATION
C     ON THE COLORING. IN DETAIL:
C
C        ICG(I) = 1: I IS C-POINT;
C        ICG(I) = 0: I IS FORCED F-POINT, I.E., IT IS AN F-POINT WHICH
C                    HAS NO STRONG CONNECTION AT ALL;
C        ICG(I) =-1: I IS "REGULAR" F-POINT, I.E., I HAS AT LEAST ONE
C                    STRONG CONNECTION TO A C-POINT;
C        ICG(I) =-2: I IS "EXCEPTIONAL" F-POINT, I.E., I HAS NO STRONG
C                    CONNECTION TO A C-POINT. FURTHERMORE, NO REGULAR
C                    F-POINT IS STRONGLY CONNECTED TO I AND ALL POINTS
C                    WITH ICG=-2 HAVE NO STRONG CONNECTION AMONG EACH
C                    OTHER. (NOTE: THESE POINTS DON'T CONTRIBUTE FROM
C                    ANY C-POINT. ALSO, IT DOES NOT MAKE SENSE TO MAKE
C                    THESE POINTS C-POINTS AS NO F-POINT CONTRIBUTES
C                    FROM THEM.)
C
C     ============ COMMENTS ON INPUT ===================================
C
C     - JA(IA(I)) ----- I=IMIN(K),...,IMAX(K)
C
C     USED AS DEFINED IN RWSRT. NOT CHANGED.
C
C     - JTR(J)--------- J=1,IW(IMAX(K)+IWS+1)-1
C     - IW(I) --------- I=IMIN(K)+IWS,...,IMAX(K)+IWS
C
C     (IWS=0, IF K=1; IWS=IMAXW(K-1)+2-IMIN(K) OTHERWISE)
C
C     USED AS DEFINED IN RWSRT. NOT CHANGED.
C
C     ============== COMMENTS ON WORK SPACE USED =======================
C
C     - ICG(I) -------- I=IMIN(K),...,IMAX(K)
C
C     USED TO DISTINGUISH F-, C- AND U- (UNDECIDED) POINTS.
C     ON ENTRY, ICG IS DEFINED TO BE 0 FOR FORCED F-POINTS AND -2
C     FOR U-PNTS. ON EXIT: SEE ABOVE.
C
C     - IFG(I) -------- I=IMIN(K),...,IMAX(K)
C
C     DEFINED TO BE A MEASURE FOR IMPORTANCE OF MAKING THE U-PNT I
C     A C-POINT. THIS MEASURE IS A VALUE BETWEEN JVAL0 AND JVALX WITH
C
C       JVAL0=IMAX(K)+NPTS+2,   NPTS=IMAX(K)-IMIN(K)+1;
C       JVALX=JVAL0+4*NTRAV+1,  NTRAV=AVVERAGE NUMBER OF STRONG TRANS-
C                                     POSE CONNECTIONS
C
C     THE HIGHER THIS VALUE, THE MORE IMPORTANT IS IT TO MAKE U-POINT I
C     A C-POINT. ALSO, THE VALUE JV=IFG(I) IS JUST THE "ORIGIN" OF
C     A LIST WHICH CONNECTS ALL POINTS HAVING THE SAME MEASURE JV OF
C     IMPORTANCE.
C
C     - ICG(II) -------- II=IMAX(K)+2,...,JVAL0-1
C     - IFG(II) -------- II=IMAX(K)+2,...,JVAL0-1
C
C     LEFT AND RIGHT STACK POINTERS IN A DOUBLY LINKED LIST. THERE
C     ARE SEVERAL SUCH LISTS, EACH OF THEM LINKING U-POINTS WHICH
C     HAVE THE SAME MEASURE OF IMPORTANCE JV TO BE MADE C-POINTS.
C     EACH OF THE VALUES JV MAY BE REGARDED AS THE "ORIGIN" OF SUCH
C     A LIST DESCRIBED IN THE FOLLOWING.
C
C     - ICG(JV) -------- JV=JVAL0,...,JVALX
C     - IFG(JV) -------- JV=JVAL0,...,JVALX
C
C     USED AS POINTERS INTO THE LISTS TO INDICATE ITS BEGINNING AND ITS
C     END: IFG(JV) POINTS TO THE FIRST POINT IN THE LIST AND ICG(JV)
C     POINTS TO THE LAST ONE. THE ROUGH PICTURE OF THE LIST WHICH
C     CORRESPONDS TO A VALUE JVAL0 <= JV <= JVALX LOOKS AS  FOLLOWS:
C
C
C                                  IFG
C       ----------------------------------------------------------
C       |                                                        |
C       |            IFG           IFG           IFG             |
C       ---> ------ ----> ------- ----> ------- ----> ------- ----
C            |    |       |     |       |     |       |     |
C            | JV |       | II1 |       | II2 |       | II3 |
C            |    |  ICG  |     |  ICG  |     |  ICG  |     |
C       ---- ------ <---- ------- <---- ------- <---- ------- <---
C       |                                                        |
C       |                          ICG                           |
C       ----------------------------------------------------------
C
C
C     HERE, II1, II2,... LOGICALLY REPRESENT THE "PHYSICAL" POINTS
C     I1, I2,... (WHICH ARE VALUES BETWEEN IMIN(K) AND IMAX(K)).
C     THE INTERCONNECTION BETWEEN THESE TWO REPRESENTATIONS OF THE
C     SAME POINTS IS GIVEN BY THE RELATION
C
C                          II := I+NPTS+1.
C
C     AN EMPTY LIST IS CHARACTERIZED BY ICG(JV)=JV, IFG(JV)=JV.
C     OBVIOUSLY, IT IS QUITE EASY TO REMOVE OR ADD POINTS TO THE LIST.
C     SUCH RE-ARRANGEMENTS ARE NECESSARY AS (IN THE COLORING PART OF
C     THE COARSENING ALGORITHM) UNDECIDED POINTS BECOME C- OR F-POINTS
C     (THEY HAVE TO BE REMOVED FROM THEIR LIST) AND THE MEASURE VALUES
C     JV OF SOME POINTS CHANGE DURING THE ALGORITHM. SUCH POINTS HAVE
C     TO BE MOVED FROM ONE LIST TO ANOTHER ONE. IN ORDER TO KEEP TRACK
C     OF THOSE POINTS WHICH HAVE TO BE RE-ARRANGED, A RESET-STACK IS
C     USED WHICH CONTAINS ALL THESE POINTS I (I.E. THEIR "LOGICAL"
C     NUMBERS II). THE POINTER IR IS USED TO POINT FROM ONE POINT IN
C     THE STACK TO THE PREVIOUS ONE.
C
C     - IR(J) --------- J=1,...,NPTS
C
C     USED BELOW AS STACK-POINTER FOR POINTS TO BE RESET IN PERFORMING
C     THE COLORING ALGORITHM. ITOP IS THE TOP-OF-STACK POINTER.
C     THE STACK IS EMPTY IF ITOP=-1.
C
C     THE GLOBAL PICTURE OF THE POINTERS IS SKETCHED IN THE FOLLOWING
C
C
C                                          LIST END
C                                       <-------------|
C                                         LIST START  |
C                                       <------------||
C                                                    ||
C                  |-- II=I+NPTS+1 ->|            IFG||ICG
C                  |                 |               ||
C                  |                 |               ||
C       IMIN(K)          IMAX(K)           JVAL0             JVALX
C          |       I        |        II      |       JV        |
C          |=======*========|========*=======|========*========|
C          |                |                |                 |
C               GRID K          LOGICAL NUM      LIST ORIGINS
C                                              (MEASURE VALUES)
C                 ||                                  |
C                 ||         IFG                      |
C                 ||--------------------------------->|
C                 |
C                 |     |---> 1 (C)
C                 |     |
C                 | ICG |---> 0 (FF)
C                 |-----|
C                       |--->-1 (F)
C                       |
C                       |--->-2 (U)
C
C
C          1                NPTS
C          |        J        |
C          |========*========|
C          |                 |
C           SHIFTED GRID K
C                   |
C                   | IR
C                   |-----> RESET STACK (ITOP = TOP-OF-STACK)
C
C
      REAL TOLD,TNEW,TIME(*)
      INTEGER IA(*),JA(*),IW(*),JTR(*),IR(*),IFG(*),ICG(*),IMIN(*),
     +        IMAX(*),IMAXW(*)
      EXTERNAL CTIME
      SAVE
C
C===> PREPARATION
C
      CALL CTIME(TOLD)
      ILO = IMIN(K)
      IHI = IMAX(K)
      IF (K.NE.1) THEN
        IWS = IMAXW(K-1)+2-ILO
      ELSE
        IWS = 0
      ENDIF
      ILO1 = ILO-1
      NPTS = IMAX(K)-IMIN(K)+1
      NPTS1 = NPTS+1
      NTRLIM = 2*(IW(IHI+IWS+1)-IW(ILO+IWS))/NPTS
      NSCNMX = 0
      JVAL0 = IHI+NPTS1+1
      JVALX = JVAL0+2*NTRLIM+1
      MDICG = MAX(MDICG,JVALX+IIAS)
      MDIR  = MAX(MDIR,NPTS)
      IF (JVALX.GT.NDICG) GOTO 9901
      IF (NPTS.GT.NDIR) GOTO 9902
C
C===> PUT INITIAL "MEASURE" FOR EACH POINT I INTO IFG(I).
C
      DO 2 I=ILO,IHI
        NSCN = IW(I+IWS+1)-IW(I+IWS)
        IF (NSCN.LE.NTRLIM) THEN
          IFG(I) = JVAL0+NSCN
        ELSE
          IFG(I) = JVALX
        ENDIF
    2 CONTINUE
C
C===> SET CIRCULARLY LINKED LISTS AND RESET-STACK TO EMPTY
C
      ITOP = -1
      DO 10 J=JVAL0,JVALX
        ICG(J) = J
        IFG(J) = J
   10 CONTINUE
C
C===> PUT ALL U-POINTS OF GRID K INTO LISTS (I.E., NO FORCED F-POINTS).
C     ADD POINTS ALWAYS TO THE END OF THEIR CORRESPONDING LIST.
C     IN THE FOLLOWING, JCNBHI DENOTES THE ACTUAL HIGHEST MEASURE VALUE
C     (AMONG U-POINTS).
C
      JCNBHI = 0
      DO 20 I=ILO,IHI
        IF(JA(IA(I)).GT.IA(I)) GOTO 15
        ICG(I) = 0
        GOTO 20
  15    ICG(I) = -2
        II = I+NPTS1
        IR(I-ILO1) = 0
        JV = IFG(I)
        IF (JV.GT.JCNBHI) JCNBHI = JV
        ICG(II) = ICG(JV)
        IFG(II) = JV
        ICG(JV) = II
        IFG(ICG(II)) = II
   20 CONTINUE
C
C     ****************
C     * PRE-COLORING *
C     ****************
C
C     PICK A U-POINT IC WITH MAXIMAL MEASURE AS GIVEN BY JCNBHI.
C     (TAKE THE FIRST FROM CORRESPONDING LIST: FIRST-IN/FIRST-OUT)
C     MAKE THAT POINT A C-POINT AND REMOVE IT FROM LISTS. THEN MAKE
C     ALL STRONG TRANSPOSE U-CONNECTIONS F-POINTS AND REMOVE THEM FROM
C     THEIR LISTS. UPDATE THE MEASURE OF IMPORTANCE FOR U-POINTS TO
C     BECOME C-POINTS AND ADD THESE POINTS TO THE RESET-STACK. FINALLY,
C     RE-ARRANGE THE LISTS BY SWEEPING THROUGH THE RESET STACK. THEN
C     PICK ANOTHER U-POINT IC.
C
C     IF LIST CORRESPONDING TO JCNBHI IS EMPTY, GO TO NEXT LOWER VALUE.
C     PRE-COLOURING IS FINISHED IF JCNBHI<=JVAL0. ALL U-POINTS LEFT AT
C     THAT TIME, WILL BE REGARDED AS F-POINTS LATER.
C
   30 IF(JCNBHI.LE.JVAL0) GOTO 100
      IIC = IFG(JCNBHI)
      IF (IIC.NE.JCNBHI) GOTO 40
      JCNBHI = JCNBHI-1
      GOTO 30
C
C===> CREATE C-POINT
C
   40 IC = IIC-NPTS1
      ICG(IC) = 1
      ICG(IFG(IIC)) = ICG(IIC)
      IFG(ICG(IIC)) = IFG(IIC)
C
C===> FOR POINT IC WITH ECCESSIVE NUMBER OF STRONG TRANSPOSE
C     CONNECTIONS: MAKE IT A C-POINT BUT LET ITS CONNECTED POINTS
C     REMAIN UNDECIDED.
C
      IF (JCNBHI.EQ.JVALX) GOTO 78
C
C===> CREATE F-POINTS AROUND ABOVE C-POINT
C
      DO 77 J=IW(IC+IWS),IW(IC+IWS+1)-1
        I = JTR(J)
        IF (ICG(I).NE.-2) GOTO 77
        ICG(I) = -1
        II = I+NPTS1
        ICG(IFG(II)) = ICG(II)
        IFG(ICG(II)) = IFG(II)
C
C===>   INCREMENT MEASURE FOR ALL STRONG U-CONNECTIONS III OF I
C       (IF NOT YET MEASURE=JVALX) AND PUT THEM ON RESET STACK (IF NOT
C       YET THERE)
C
        DO 76 JJ=IA(I)+1,JA(IA(I))
          III = JA(JJ)
          IF (ICG(III).NE.-2.OR.IFG(III).GE.JVALX) GOTO 76
          IFG(III) = IFG(III)+1
          IIII = III-ILO1
          IF (IR(IIII).NE.0) GOTO 76
          IR(IIII) = ITOP
          ITOP = IIII
   76   CONTINUE
   77 CONTINUE
C
C===> DECREMENT MEASURE FOR ALL STRONG U-CONNECTIONS I OF IC
C     AND PUT THEM ON RESET-STACK (IF NOT YET THERE)
C
   78 DO 87 J=IA(IC)+1,JA(IA(IC))
        I = JA(J)
        IF (ICG(I).NE.-2) GOTO 87
        IFG(I) = IFG(I)-1
        II = I-ILO1
        IF (IR(II).NE.0) GOTO 87
        IR(II) = ITOP
        ITOP = II
   87 CONTINUE
C
C===> REARRANGE THE LISTS BY SWEEPING THROUGH RESET-STACK.
C     THEN GO BACK TO PICK ANOTHER U-POINT IC.
C
      IN = ITOP
      ITOP = -1
C
   90 IF (IN.LE.0) GOTO 30
      I = IN+ILO1
      II = I+NPTS1
      IF (ICG(I).NE.-2) GOTO 95
      IFG(ICG(II)) = IFG(II)
      ICG(IFG(II)) = ICG(II)
      JV = IFG(I)
      IF (JV.GT.JCNBHI) JCNBHI = JV
      ICG(II) = ICG(JV)
      IFG(II) = JV
      ICG(JV) = II
      IFG(ICG(II)) = II
   95 IIP = IN
      IN = IR(IN)
      IR(IIP) = 0
      GOTO 90
C
  100 CALL CTIME(TNEW)
      TIME(2) = TIME(2)+TNEW-TOLD
      RETURN
C
C===> ERROR MESSAGES
C
 9901 WRITE (IUM,9910)
      IERR = 6
      RETURN
C
 9902 WRITE (IUM,9920)
      IERR = 4
      RETURN
C
 9910 FORMAT (' *** ERROR IN PCOL: NDIG TOO SMALL ***')
 9920 FORMAT (' *** ERROR IN PCOL: NDU TOO SMALL ***')
      END
C
C.......................................................................
C
C     PWINT                                             SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE PWINT (K,EWT2,ICHK,MMAX,A,IA,JA,IW,IMIN,IMAX,IMINW,
     +                  IMAXW,IFG,ICG,NCOLOR,TIME,IERR,IROW0,NCOLX,
     +                  NDA,NDJA,NDICG,NDIA,IUM,MDA,MDJA,IAJAS)
C
C     SET UP FINAL COARSER GRID K+1 AND INTERPOLATION FORMULA FROM
C     GRID K+1 TO GRID K. THIS IS THE VERSION AS DESCRIBED IN RUGE/
C     STUEBEN (BRISTOL). PWINT ASSUMES THE GRID TO BE PRE-COLORED
C     BY SUBROUTINE PCOL.
C
C     ON EXIT, A, JA AND IW ARE SET TO CONTAIN THE INTERPOLATION
C     WEIGHTS AND CORRESPONDING POINTERS AS REQUIRED IN THE SOLUTION
C     PHASE OF AMG1R5. ALSO, ICG(I) (IMIN(K)<=I<=IMAX(K)) ARE SET TO
C     THEIR FINAL VALUES, EXCEPT FOR THOSE I WITH ICG(I)<0.
C
C     ============ COMMENTS ON INPUT ===================================
C
C     - JA(IA(I)) ----- I=IMIN(K),...,IMAX(K)
C
C     ASSUMED TO POINT TO THE LAST STRONG CONNECTION OF POINT I (OR TO
C     IA(I) IF THERE IS NO SUCH CONNECTION). ON EXIT, RESET TO ORIGINAL
C     VALUES (I.E. JA(IA(I))=I).
C
C     - ICG(I) -------- I=IMIN(K),...,IMAX(K)
C
C     ASSUMED TO CONTAIN INFORMATION ON PRE-COLORING:
C
C       ICG(I)>0: I IS C-POINT
C       ICG(I)=0: I IS FORCED F-POINT (I.E. I HAS NO STRONG CONNNECTION)
C       ICG(I)<0: I IS F-POINT WITH (AT LEAST) ONE STRONG CONNECTION.
C
C     ============== COMMENTS ON WORK SPACE USED =======================
C
C     - IFG(I) -------- I=IMIN(K),...,IMAX(K)
C
C     IS USED FOR SEVERAL PURPOSES. IN PARTICULAR, TO DISTINGUISH
C     INTERPOLATORY AND NON-INTERPOLATORY POINTS.
C
C     - NCOLOR(I) ----- I=1,...,#POINTS ON GRID K
C
C     IS SET TO F-POINT-COLORS TO BE USED LATER IN SUBROUTINE OPDFN.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION A(*)
      REAL TOLD,TNEW,TIME(*)
      INTEGER IA(*),JA(*),IW(*),IMIN(*),IMAX(*),IMINW(*),IMAXW(*),
     +        IFG(*),ICG(*),NCOLOR(*)
      EXTERNAL CTIME
      SAVE
C
      CALL CTIME(TOLD)
      NDAJA = MIN(NDA,NDJA)
      NCOUNT = 0
      IF (K.EQ.1) THEN
        IMINW(1) = 1
        IW(1) = IA(IMAX(1)+1)
      ENDIF
C
      ILO = IMIN(K)
      IHI = IMAX(K)
      DO 12 I=ILO,IHI
        IFG(I) = 0
   12 CONTINUE
C
C===> SWEEP OVER F-POINTS I WHICH HAVE AT LEAST ONE STRONG CONNECTION
C
      IBLCK = IMINW(K)
      DO 400 I=ILO,IHI
        IF (ICG(I).GE.0) GOTO 400
        JLO = IA(I)+1
        JHI = JA(IA(I))
        EWT2I = EWT2/A(JLO)
        NCONDC = 0
C
C===>   INITIALIZE "BLOCK" OF INTERPOLATION WEIGHTS OF POINT I
C
   30   JW0 = IW(IBLCK)
        JWX = JW0
        IF (JWX.GT.NDAJA) GOTO 2000
        DO 20 J=JLO,JHI
          II = JA(J)
          IF (ICG(II).LE.0) GOTO 20
          A(JWX)  = A(J)
          JA(JWX) = II
          IFG(II) = JWX
          JWX = JWX+1
          IF (JWX.GT.NDAJA) GOTO 2000
   20   CONTINUE
        A(JWX)  = A(IA(I))
        JA(JWX) = I
        MDA  = MAX(MDA,JWX+IAJAS)
        MDJA = MAX(MDJA,JWX+IAJAS)
        IFG(I) = JWX
C
C===>   SWEEP OVER STRONGLY CONNECTED F-POINTS. THESE MUST BE "COVERED"
C       BY A TOTAL WEIGHT DEFINED BY EWT2. IF AN F-POINT HAS NO STRONG
C       CONNECTIONS, REGARD IT TO BE COVERED, BUT DO NOT DISTRIBUTE
C       THE CORRESPONDING WEIGHT (ERROR AT SUCH A POINT CAN BE ASSUMED
C       TO BE VERY SMALL!).
C
        DO 150 J=JLO,JHI
          II = JA(J)
          IF (ICG(II).GE.0) GOTO 150
C
C===>     COMPUTE DEPENDENCE ON SET OF INTERPOLATION POINTS
C
          S  = 0.0D0
          SI = 0.0D0
          JJLO = IA(II)+1
          JJHI = IA(II+1)-1
          DO 110 JJ=JJLO,JJHI
            IF (IFG(JA(JJ)).LT.JW0) GOTO 110
            IF (JA(JJ).EQ.I) SI = A(JJ)
            S = S+A(JJ)
  110     CONTINUE
          IF (ICHK.EQ.2) GOTO 111
          IF (S.EQ.0.0D0) THEN
            A(JWX) = A(JWX)+A(J)
            GOTO 150
          ELSE
            GOTO 135
          ENDIF
C
C===>     CHECK DEPENDENCE ON SET OF INTERPOLATION POINTS
C
  111     IF (S-SI.LE.EWT2I*A(J)*A(JJLO)) GOTO 135
C
C===>     DEPENDENCE TOO SMALL: IF THERE IS NOT YET A CONDITIONAL
C         C-POINT, MAKE II SUCH A POINT AND RESTART THE
C         PROCESS FOR DEFINING INTERPOLATION WEIGHTS FOR POINT I.
C         OTHERWISE MAKE I ITSELF A C-POINT AND LEAVE II AN F-POINT.
C
          IF (NCONDC.EQ.0) THEN
            NCOUNT = NCOUNT+1
            NCONDC = 1
            IP = II
            ICGP = ICG(II)
            ICG(II) = 1
            GOTO 30
          ELSE
            ICG(I) = 1
            ICG(IP) = ICGP
            DO 120 JJ=JW0,JWX
              IFG(JA(JJ)) = 0
  120       CONTINUE
            GOTO 400
          ENDIF
C
C===>     DISTRIBUTE THE WEIGHT OF POINT II
C
  135     WW = A(J)/S
          DO 140 JJ=JJLO,JJHI
            IF (IFG(JA(JJ)).GE.JW0)
     +          A(IFG(JA(JJ))) = A(IFG(JA(JJ)))+A(JJ)*WW
  140     CONTINUE
  150   CONTINUE
C
C===>   ALL NECESSARY POINTS ARE COVERED. NOW DISTRIBUTE WEIGHTS FROM
C       WEAK CONNECTIONS OF POINT I (ANALOGOUS AS ABOVE)
C
        DO 190 J=JHI+1,IA(I+1)-1
          II = JA(J)
          IF (ICG(II).EQ.0) GOTO 190
          S = 0.0D0
          JJLO = IA(II)+1
          JJHI = IA(II+1)-1
          DO 160 JJ=JJLO,JJHI
            IF (IFG(JA(JJ)).GE.JW0) S = S+A(JJ)
  160     CONTINUE
          IF (S.EQ.0.0D0) THEN
            A(JWX) = A(JWX)+A(J)
          ELSE
            WW = A(J)/S
            DO 170 JJ=JJLO,JJHI
              IF (IFG(JA(JJ)).GE.JW0)
     *           A(IFG(JA(JJ))) = A(IFG(JA(JJ)))+A(JJ)*WW
  170       CONTINUE
          ENDIF
  190   CONTINUE
C
        ICG(I) = -IBLCK
        IBLCK = IBLCK+1
        IW(IBLCK) = JWX+1
  400 CONTINUE
C
C===> SET ICG; RESET JA(IA(I)); CHECK SIZE OF COARSEST GRID
C
      IC = IHI
      IMIN(K+1) = IC+1
      DO 900 I=ILO,IHI
        JA(IA(I)) = I
        IF (ICG(I).LE.0) GOTO 900
        IC = IC+1
        ICG(I) = IC
        IF (IC.LT.NDICG) ICG(IC) = 0
  900 CONTINUE
      IMAX(K+1) = IC
C
      NPTS = IHI-ILO+1
      NPTSC = IMAX(K+1)-IMIN(K+1)+1
      IF (NPTSC.EQ.1)                  MMAX = K+1
      IF (NPTSC.EQ.1.AND.IROW0.LT.2)   MMAX = K
      IF (NPTSC.EQ.NPTS.OR.NPTSC.EQ.0) MMAX = K
      IF (K.GE.MMAX) GOTO 1000
      IF (IC.GE.NDIA) GOTO 1700
      IF (IC.GE.NDICG) GOTO 1800
C
C===> RE-ARRANGE A
C
 1000 IBLCK1 = IMINW(K)
      JWPOS = IW(IBLCK1)
      DO 950 I=ILO,IHI
        IF (ICG(I).GE.0) GOTO 950
        IBLCK = -ICG(I)
        J1 = IW(IBLCK)
        J2 = IW(IBLCK+1)-1
        IF (J2.LE.J1) THEN
          ICG(I) = 0
        ELSE
          ICG(I) = -IBLCK1
          IW(IBLCK1) = JWPOS
          SCALE = -1.0D0/A(J2)
          DO 920 J=J1,J2-1
            A(JWPOS)  = A(J)*SCALE
            JA(JWPOS) = ICG(JA(J))
            JWPOS = JWPOS+1
  920     CONTINUE
          IBLCK1 = IBLCK1+1
        ENDIF
  950 CONTINUE
      IMAXW(K) = IBLCK1-1
      IW(IBLCK1) = JWPOS
C
C===> STORE TYPE OF POINTS (I.E. C, F OR FF)  ON VECTOR NCOLOR
C
      IS = 1-ILO
      DO 960 I=ILO,IHI
        NC = ICG(I)
        IF (NC.LT.0) THEN
          NCOLOR(I+IS) = 1
        ELSEIF (NC.GT.0) THEN
          NCOLOR(I+IS) = 2
        ELSE
          NCOLOR(I+IS) = 3
        ENDIF
  960 CONTINUE
      NCOLX = 1
C
C===> EXIT
C
      WRITE (IUM,9000) K,NCOUNT
      CALL CTIME(TNEW)
      TIME(4) = TIME(4)+TNEW-TOLD
      RETURN
C
C===> ERROR MESSAGES
C
 1700 WRITE (IUM,9030)
      IERR = 2
      RETURN
C
 1800 WRITE (IUM,9040)
      IERR = 6
      RETURN
C
 2000 IF (NDA.LE.NDJA) THEN
        WRITE (IUM,9020)
        IERR = 1
      ELSE
        WRITE (IUM,9050)
        IERR = 3
      ENDIF
      RETURN
C
 9000 FORMAT (' INTERPOLATION OPERATOR NO.',I3,' COMPLETED. C-POINTS',
     +        ' ADDED IN PWINT:',I4)
 9020 FORMAT (' *** ERROR IN PWINT: NDA TOO SMALL ***')
 9030 FORMAT (' *** ERROR IN PWINT: NDIA TOO SMALL ***')
 9040 FORMAT (' *** ERROR IN PWINT: NDIG TOO SMALL ***')
 9050 FORMAT (' *** ERROR IN PWINT: NDJA TOO SMALL ***')
      END
C
C.......................................................................
C
C     OPDFN                                             SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE OPDFN(K,IERR,MMAX,A,IA,JA,IW,IMIN,IMAX,IMINW,IMAXW,
     +                 ICG,IFG,NCOLOR,NSTCOL,NCOLX,TIME,NDA,NDJA,
     +                 IUM,MDA,MDJA,IAJAS)
C
C     THIS SUBROUTINE CONSTRUCTS THE CG-OPERATOR A(K) (K>1).
C
C     ONE ROW IS CONSTRUCTED AT A TIME, AND THE ROW STRUCTURE
C     (I.E., WHICH CONNECTION GOES WHERE) MUST BE DETERMINED
C     AT THE SAME TIME.  IN ORDER TO AVOID SEARCHES THROUGH
C     THE CURRENT ROW TO DETERMINE IF A POSITION FOR A CONNECTION
C     HAS ALREADY BEEN DEFINED, ICG (FOR LEVEL K) AND IFG (FOR LEVEL
C     K-1) ARE USED AUXILIARILY. FURTHERMORE, TO SPEED UP COMPUTATION,
C     THE TRANSPOSE OF INTERPOLATION IS TEMPORARILY STORED ON A (AT
C     THE END OF THE AVAILABLE SPACE). IFG (FOR LEVEL K) IS USED
C     AS POINTER TO THE CORRESPONDING ROWS. MORE DETAILS: SEE BELOW.
C
C     ============ COMMENTS ON INPUT ===================================
C
C     - IMIN(K), IMAX(K)
C
C     FIRST/LAST NUMBER OF GRID POINTS ON GRID K.
C
C     - IMINW(K-1), IMAXW(K-1)
C
C     FIRST/LAST "BLOCK" NUMBER OF INTERPOLATION WEIGHTS CONTRIBUTING
C     IN INTERPOLATION TO GRID K-1.
C
C     - IW(IBL) ----- IBL=IMINW(K-1),IMAXW(K-1)+1
C     - A(J) -------- J=JLO,JHI      (JLO=IW(IMINW(K-1)))
C     - JA(J) ------- J=JLO,JHI      (JHI=IW(IMAXW(K-1)+1)-1)
C
C     ARE ASSUMED TO CONTAIN THE WEIGHTS OF INTERPOLATION ALONG WITH
C     THE NECESSARY POINTERS. IW(IBL) POINTS TO THE FIRST ENTRY OF
C     "BLOCK" NUMBER IBL IN A (CF. BELOW).
C
C     - ICG(IF) ----- IF=IMIN(K-1),IMAX(K-1)
C
C     ARE ASSUMED TO BE SET TO THEIR FINAL VALUES:
C
C       ICG(IF)>0: IF IS C-POINT OF GRID K-1 AND ICG(IF) POINTS JUST
C                  TO THE CORRESPONDING POINT ON GRID K;
C       ICG(IF)=0: IF IS F-POINT OF GRID K-1 WITHOUT ANY CONTRIBUTION
C                  IN INTERPOLATION FROM GRID K;
C       ICG(IF)<0: IF IS F-POINT OF GRID K-1; IBL=-ICG(IF) POINTS JUST
C                  TO THE "BLOCK" OF INTERPOLATION WEIGHTS. THAT MEANS,
C                  JW(J) (IW(IBL)<=J<=IW(IBL+1)-1) POINTS TO THE POINTS
C                  (ON GRID K) WHICH CONTRIBUTE IN INTERPOLATION TO IF.
C                  THE CORRESPONDING WEIGHTS ARE STORED IN A(J).
C
C     - IFG(IF) ----- IF=IMIN(K-1),IMAX(K-1)
C
C     WORK SPACE (SEE BELOW). IFG(IF) IS ASSUMED TO BE > -IMIN(K).
C
C     - IFG(IC) ----- IC=IMIN(K),IMAX(K)+1
C
C     WORK SPACE (SEE BELOW). THE CONTENTS OF IFG(IC) IS ARBITRARY.
C
C     - ICG(IC) ----- IC=IMIN(K),IMAX(K)
C
C     WORK SPACE (SEE BELOW). ICG(IC) IS ASSUMED TO BE ZERO.
C
C     ============== COMMENTS ON WORK SPACE USED =======================
C
C     - A(J) -------- J=JJLO,JJHI
C     - JA(J) ------- J=JJLO,JJHI
C
C     (JJLO=MIN(NDA,NDJA)-IW(IMAXW(K-1)+1)+IW(IMINW(K-1)+1),
C      JJHI=MIN(NDA,NDJA), I.E. THE SPACE OF LENGTH OF INTERPOLATION
C      (K-1) AT THE END OF A AND JA, RESPECTIVELY.)
C
C     IS USED TO STORE THE TRANSPOSE OF INTERPOLATION W(K-1). THIS IS
C     DONE IN ORDER TO SPEED UP THE COMPUTATION OF THE OPERATOR A(K).
C     IF DURING ASSEMBLAGE OF A(K) THIS WORK SPACE IS REQUIRED BY A(K),
C     PROCESSING CONTINUES USING THE NOT YET REDEFINED PART OF THE WORK
C     SPACE AND RECALCULATING THE LOST ENTRIES EACH TIME THEY ARE
C     NEEDED. THUS THE CALCULATION SLOWS DOWN.
C
C     - IFG(IC) ----- IC=IMIN(K),IMAX(K)+1
C
C     IS USED AS POINTER FOR THE TRANSPOSE OF INTERPOLATION: THE COARSE-
C     GRID POINT IC CONTRIBUTES IN INTERPOLATION TO THE FINE-GRID POINTS
C     IF=JA(J) (IFG(IC)<=J<=IFG(IC+1)-1) WITH WEIGHT A(J). THE
C     CONTRIBUTION TO ITSELF (BY THE WEIGHT 1.0) IS NOT CONTAINED.
C
C     - ICG(IC) ----- IC=IMIN(K),IMAX(K)
C
C     IS USED FOR SEVERAL PURPOSES. IN ASSEMBLING THE CG-OPERATOR A(K),
C     IT SERVES AS POINTER TO POSITIONS IN A(K) WHICH HAVE ALREADY BEEN
C     DEFINED: IF THE CURRENT ROW CORRESPONDS TO POINT IC1, AND A
C     CONNECTION TO ANOTHER CG-POINT IC2 HAS JUST BEEN FOUND, THEN IF
C     ICG(IC2)<IA(IC1), THE CORRESPONDING ENTRY IN ROW IC1 OF A(K)
C     HAS NOT YET BEEN DEFINED. OTHERWISE, ICG(IC2) POINTS TO THE
C     LOCATION FOR THAT ENTRY. (ALSO SEE IFG BELOW.)
C
C     - IFG(IF) ----- IF=IMIN(K-1),IMAX(K-1)
C
C     IN ASSEMBLING THE CG-OPERATOR A(K), THIS VECTOR CONTAINS INFORMAT-
C     ION ON WHETHER THE EXISTENCE OF ENTRIES IN A(K) HAS TO BE CHECKED
C     OR NOT.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION A(*)
      REAL TOLD,TNEW,TIME(*)
      INTEGER IA(*),JA(*),IW(*),IMIN(*),IMAX(*),IMINW(*),IMAXW(*),
     +        ICG(*),IFG(*),NCOLOR(*),NSTCOL(*)
      EXTERNAL CTIME
      SAVE
C
C===> EXTEND A, JA TO STORE TRANSPOSE OF INTERPOLATION
C
      CALL CTIME(TOLD)
      NDAJA = MIN(NDA,NDJA)
      IMINW(K) = IMAXW(K-1)+1
      IF (K.GT.MMAX) GOTO 8000
      JPOS = IW(IMINW(K))
      DO 200 J=IW(IMINW(K-1)),IW(IMAXW(K-1)+1)-1
        ICG(JA(J)) = ICG(JA(J))+1
  200 CONTINUE
C
      MDA  = MAX(MDA,JPOS+JPOS-IW(IMINW(K-1))+IAJAS)
      MDJA = MAX(MDJA,JPOS+JPOS-IW(IMINW(K-1))+IAJAS)
      IFG(IMIN(K)) = NDAJA-JPOS+IW(IMINW(K-1))+1
      IF (IFG(IMIN(K)).LE.JPOS) GOTO 9900
      DO 220 IC=IMIN(K),IMAX(K)
        IFG(IC+1) = IFG(IC)+ICG(IC)
        ICG(IC) = IFG(IC)
  220 CONTINUE
C
      DO 250 IF=IMIN(K-1),IMAX(K-1)
        IF (ICG(IF).GE.0) GOTO 250
        IBL = -ICG(IF)
        DO 260 J=IW(IBL),IW(IBL+1)-1
          IC = JA(J)
          A(ICG(IC))=A(J)
          JA(ICG(IC))=IF
          ICG(IC)=ICG(IC)+1
  260   CONTINUE
  250 CONTINUE
C
      DO 300 IC=IMIN(K),IMAX(K)
        ICG(IC) = 0
  300 CONTINUE
C
C===> SWEEP OVER ALL CG-POINTS IC TO ASSEMBLE ROWS OF CG MATRIX
C
      ISTTI = 1
      DO 100 IF=IMIN(K-1),IMAX(K-1)
        IF (ICG(IF).LE.0) GOTO 100
        IC = ICG(IF)
        IF (JPOS.GT.NDAJA) GOTO 9900
        IALOW = JPOS
        ICG(IC) = JPOS
        A(JPOS)  = A(IA(IF))
        JA(JPOS) = IC
        JPOS = JPOS+1
C
C       ----------------------------------------------
C       | SEARCH FOR C-C-C-C AND C-C-F-C CONNECTIONS |
C       ----------------------------------------------
C
        DO 25 JF1=IA(IF)+1,IA(IF+1)-1
          IF1 = JA(JF1)
          IC1 = ICG(IF1)
          IF (IC1) 11,25,20
C
C===>     IF1 IS F-POINT: SWEEP OVER C-C-F-C CONNECTIONS
C
   11     IFG(IF1) = -IC
          DO 15 JF2=IW(-IC1),IW(-IC1+1)-1
            IC2 = JA(JF2)
            IF (ICG(IC2).GE.IALOW) GOTO 10
            IF (JPOS.GT.NDAJA) GOTO 9900
            ICG(IC2) = JPOS
            A(JPOS)  = A(JF1)*A(JF2)
            JA(JPOS) = IC2
            JPOS = JPOS+1
            GOTO 15
   10       A(ICG(IC2)) = A(ICG(IC2))+A(JF1)*A(JF2)
   15     CONTINUE
          GOTO 25
C
C===>     IF1 IS C-POINT: C-C-C-C CONNECTION
C
   20     IF (ICG(IC1).GE.IALOW) GOTO 23
          IF (JPOS.GT.NDAJA) GOTO 9900
          ICG(IC1) = JPOS
          A(JPOS)  = A(JF1)
          JA(JPOS) = IC1
          JPOS = JPOS+1
          GOTO 25
   23     A(ICG(IC1)) = A(ICG(IC1))+A(JF1)
   25   CONTINUE
        IST = IMIN(K-1)-1
C
C        ----------------------------------------------
C        | SEARCH FOR C-F-C-C AND C-F-F-C CONNECTIONS |
C        ----------------------------------------------
C
        DO 95 JF1=IFG(IC),IFG(IC+1)-1
          IF (JF1.GE.JPOS) THEN
            IF1 = JA(JF1)
            WJF1 = A(JF1)
            IST = IF1
          ELSE
            ISTTI = 0
            DO 120 JF2=IST+1,IMAX(K-1)
              IF (ICG(JF2).GE.0) GOTO 120
              IBL = -ICG(JF2)
              DO 110 JB=IW(IBL),IW(IBL+1)-1
                JC3 = JA(JB)
                IF(JC3.NE.IC) GOTO 110
                IF1 = JF2
                WJF1 = A(JB)
                GOTO 130
  110         CONTINUE
  120       CONTINUE
C
C===>       ERROR EXIT
C
            WRITE (IUM,9940) K-1
            IERR = 22
            RETURN
  130       IST = IF1
          ENDIF
          DO 90 JF2=IA(IF1),IA(IF1+1)-1
            IF2 = JA(JF2)
            IC2 = ICG(IF2)
            WW = WJF1*A(JF2)
            IF (IC2) 35,90,50
C
C===>       IF2 IS F-POINT: SWEEP OVER C-F-F-C CONNECTIONS
C
   35       IF (IFG(IF2).EQ.-IC) GOTO 70
            IFG(IF2) = -IC
            DO 40 JF3=IW(-IC2),IW(-IC2+1)-1
              IC3 = JA(JF3)
              IF (ICG(IC3).GE.IALOW) GOTO 30
              IF (JPOS.GT.NDAJA) GOTO 9900
              ICG(IC3) = JPOS
              A(JPOS) = WW*A(JF3)
              JA(JPOS) = IC3
              JPOS = JPOS+1
              GOTO 40
   30         A(ICG(IC3)) = A(ICG(IC3))+WW*A(JF3)
   40       CONTINUE
            GOTO 90
C
C===>       IF2 HAS BEEN ENCOUNTERED BEFORE; DO NOT CHECK POSITIONS!
C
   70       DO 80 JF3=IW(-IC2),IW(-IC2+1)-1
              IADRS = ICG(JA(JF3))
              A(IADRS) = A(IADRS)+WW*A(JF3)
   80       CONTINUE
            GOTO 90
C
C===>       IF2 IS C-POINT: C-F-C-C CONNECTION
C
   50       IF (ICG(IC2).GE.IALOW) GOTO 60
            IF (JPOS.GT.NDAJA) GOTO 9900
            ICG(IC2) = JPOS
            A(JPOS) = WW
            JA(JPOS) = IC2
            JPOS = JPOS+1
            GOTO 90
   60       A(ICG(IC2)) = A(ICG(IC2))+WW
   90     CONTINUE
   95   CONTINUE
        IA(IC+1) = JPOS
  100 CONTINUE
      MDA  = MAX(MDA,IA(IMAX(K)+1)-1+IAJAS)
      MDJA = MAX(MDJA,IA(IMAX(K)+1)-1+IAJAS)
C
C===> WARNING FOR STORAGE SHORTAGE
C
      IF (ISTTI.NE.1) THEN
        K1 = K-1
        WRITE (IUM,9930) K1
        IERR = -1
      ENDIF
      WRITE (IUM,9000) K
C
C===> SET UP LINKED LIST FOR RELAXATION ON GRID K-1
C
 8000 IA(IMIN(K)) = IW(IMINW(K))
      ILO = IMIN(K-1)
      IHI = IMAX(K-1)
      ILO1 = ILO-1
      NPTS = IHI-ILO1
      IST = 100000000
      DO 8100 ICOL=1,NCOLX
        DO 8040 I=NPTS,1,-1
          IF (NCOLOR(I).NE.ICOL) GOTO 8040
          ICG(I+ILO1) = -IST
          IST=I+ILO1
 8040   CONTINUE
 8100 CONTINUE
      NSTCOL(K-1)=IST
C
C===> EXIT / ERROR MESSAGES
C
      CALL CTIME(TNEW)
      TIME(6) = TIME(6)+TNEW-TOLD
      RETURN
C
 9900 IF (NDA.LE.NDJA) THEN
        WRITE (IUM,9910)
        IERR = 1
      ELSE
        WRITE (IUM,9920)
        IERR = 3
      ENDIF
      RETURN
C
 9000 FORMAT (' COARSE  GRID  OPERATOR NO.',I3,' COMPLETED')
 9910 FORMAT (' *** ERROR IN OPDFN: NDA TOO SMALL ***')
 9920 FORMAT (' *** ERROR IN OPDFN: NDJA TOO SMALL ***')
 9930 FORMAT (' --- WARNG: UNABLE TO STORE TRANSPOSE OF INTERPOLATION',
     +        ' ON GRID ',I2,' DURING '/
     +        '            EXECUTION OF OPDFN, BECAUSE NDA OR NDJA ',
     +        'TOO SMALL.'/
     +        '            SETUP COMPUTATION IS SLOWING DOWN.')
 9940 FORMAT (' *** ERROR IN OPDFN: INTERPOLATION ENTRY MISSING ON GRID'
     +        ,I3)
      END
C
C.......................................................................
C
C     SETIFG                                              SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE SETIFG(IMIN,IMAX,ICG,IFG,NSTCOL,LEVELS,TIME)
C
C     SET "INVERSE" POINTER IFG
C
      REAL TOLD,TNEW,TIME(*)
      INTEGER IMIN(*),IMAX(*),ICG(*),IFG(*),NSTCOL(*)
      EXTERNAL CTIME
      SAVE
C
      CALL CTIME(TOLD)
      DO 10 I=IMIN(1),IMAX(LEVELS-1)
        IF (ICG(I).GT.0) IFG(ICG(I))=I
   10 CONTINUE
      IB = 1
      DO 30 K=1,LEVELS-1
        IST = NSTCOL(K)
   20   IF (IST.GE.100000000) GOTO 30
        IFG(IB) = IST
        IST = -ICG(IST)
        IB = IB+1
        GOTO 20
   30 CONTINUE
      CALL CTIME(TNEW)
      TIME(8) = TIME(8)+TNEW-TOLD
      RETURN
      END
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     AMG1R5 SOLUTION-SUBROUTINES
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C.......................................................................
C
C     SOLVE                                                SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE SOLVE(MADAPT,NCYC,NRD,NSOLCO,NRU,IOUT,IERR,A,U,F,IA,
     +                 JA,IW,EPS,IMIN,IMAX,IMINW,IMAXW,ICG,IFG,NSTCOL,
     +                 IARR,TIME,NCYC0,IROW0,LEVELS,NDA,NDJA,
     +                 NDU,NDF,MDA,MDJA,MDU,MDF,IUP,IUM,RESI,RES0,RES)
C
C     SOLUTION PHASE OF AMG1R5
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION A(*),U(*),F(*),RESI(*)
      REAL TIME(*)
      INTEGER NRDTYP(10),NRUTYP(10),IA(*),JA(*),IW(*),IMIN(*),IMAX(*),
     +        IMINW(*),IMAXW(*),ICG(*),IFG(*),NSTCOL(*),IARR(*)
      EXTERNAL CG,CYC,IDEC,RESID,USAVE
      SAVE
C
C===> TEST OF AVAILABLE STORAGE
C
      IF (NDU.LT.IMAX(LEVELS)) THEN
        WRITE (IUM,9050)
        IERR = 4
        RETURN
      ENDIF
      IF (NDF.LT.IMAX(LEVELS)) THEN
        WRITE (IUM,9060)
        IERR = 5
        RETURN
      ENDIF
C
      M = LEVELS
      NCYC0 = 0
      DO 5 N=11,20
        TIME(N) = 0.0
5     CONTINUE
      IF (EPS.NE.0.D0) THEN
        EPSI = EPS
      ELSE
        EPSI = 1.D-12
      ENDIF
C
C===> DECOMPOSE MADAPT
C
      IF (MADAPT.NE.0) THEN
        CALL IDEC(MADAPT,2,NDIGIT,IARR)
        MSEL = IARR(1)
        IF (MSEL.EQ.2) THEN
          IF (IARR(2).NE.0) THEN
            FAC = DBLE(IARR(2))
            DO 8 I=1,100
              FAC = FAC/10.D0
              IF (FAC.LE.1.D0) GOTO 9
8           CONTINUE
          ELSE
            FAC = 0.7D0
          ENDIF
        ENDIF
      ELSE
        MSEL = 2
        FAC = 0.7D0
      ENDIF
C
C===> DECOMPOSE NCYC
C
9     IF (NCYC.NE.0) THEN
        CALL IDEC(ABS(NCYC),4,NDIGIT,IARR)
        IGAM   = SIGN(IARR(1),NCYC)
        ICGR   = IARR(2)
        ICONV  = IARR(3)
        NCYCLE = IARR(4)
        IF (NCYCLE.EQ.0) RETURN
      ELSE
        IGAM   = 1
        ICGR   = 0
        ICONV  = 1
        NCYCLE = 10
      ENDIF
C
C===> SET EPSI ACCORDING TO CONVERGENCE CRITERION GIVEN BY ICONV
C
      IF (ICONV.NE.3) THEN
        IF (ICONV.EQ.4) THEN
        AMA = 0.D0
        DO 6 I=1,IMAX(1)
          AMA = MAX(AMA,A(IA(I)))
6       CONTINUE
        EPSI = EPSI*AMA
        ENDIF
      ELSE
        FMAX = 0.D0
        DO 7 I=1,IMAX(1)
          FMAX = MAX(FMAX,ABS(F(I)))
7       CONTINUE
        EPSI = EPSI*FMAX
      ENDIF
C
C===> DECOMPOSE NRD
C
      IF (NRD.NE.0) THEN
        CALL IDEC(NRD,9,NDIGIT,NRDTYP)
        NRDX   = NRDTYP(2)
        NRDLEN = NDIGIT-2
        DO 10 I=1,NRDLEN
          NRDTYP(I) = NRDTYP(I+2)
10      CONTINUE
      ELSE
        NRDX      = 1
        NRDLEN    = 2
        NRDTYP(1) = 3
        NRDTYP(2) = 1
      ENDIF
C
C===> DECOMPOSE NRU
C
      IF (NRU.NE.0) THEN
        CALL IDEC(NRU,9,NDIGIT,NRUTYP)
        NRUX   = NRUTYP(2)
        NRULEN = NDIGIT-2
        DO 40 I=1,NRULEN
          NRUTYP(I) = NRUTYP(I+2)
40      CONTINUE
      ELSE
        NRUX      = 1
        NRULEN    = 2
        NRUTYP(1) = 3
        NRUTYP(2) = 1
      ENDIF
C
C===> DECOMPOSE NSOLCO
C
      IF (NSOLCO.NE.0) THEN
        CALL IDEC(NSOLCO,2,NDIGIT,IARR)
        NSC  = IARR(1)
        NRCX = IARR(2)
C
C===> IN CASE OF YALE-SMP COARSE GRID SOLUTION, DON'T USE COARSEST
C     GRID WITH LESS THAN 10 POINTS
C
        IF (NSC.EQ.2) THEN
          DO 50 I=M,1,-1
            L = I
            IF (IMAX(I)-IMIN(I).GE.9) GOTO 60
50        CONTINUE
60        M = I
          LEVELS = I
        ENDIF
      ELSE
        NSC  = 1
        NRCX = 0
      ENDIF
C
C===> CYCLING
C
100   IF (IOUT.NE.0) THEN
        CALL RESID(1,RES0,A,U,F,IA,JA,IW,IMIN,IMAX,IMINW)
        IF (IOUT.EQ.3) THEN
          WRITE (IUP,9005)
          WRITE (IUP,9000) RES0
        ENDIF
        RESOLD = RES0
      ENDIF
C
      DO 120 ITER=1,NCYCLE
        CALL USAVE(1,ICGR,U,IMIN,IMAX,NDU,M,TIME,IERR,IUM,MDU,NDF,MDF)
        CALL CYC(1,NRDX,NRDTYP,NRDLEN,NRCX,NRUX,NRUTYP,NRULEN,
     +           IGAM,A,U,F,IA,JA,IW,IMIN,IMAX,IMINW,IMAXW,
     +           IFG,ICG,NSTCOL,IARR,TIME,IROW0,M,IUM,IERR,ITER,
     +           NSC,NDA,NDJA,MDA,MDJA,MSEL,FAC,RESI,LEVELS)
        CALL CG(1,ICGR,ITER,A,U,F,IA,JA,IW,IMIN,IMAX,IMINW,M,TIME,
     +          IERR,IUM)
        IF (IERR.GT.0) RETURN
        IF (ITER.EQ.1) THEN
          MFIRST = M
          IF (IOUT.EQ.3) WRITE (IUP,9040) M
        ELSEIF (IOUT.EQ.3.AND.M.NE.MFIRST) THEN
          MFIRST = M
          WRITE (IUP,9040) M
        ENDIF
        IF (IOUT.EQ.3.OR.ICONV.NE.1) CALL RESID(1,RES,A,U,F,IA,JA,IW,
     +                                          IMIN,IMAX,IMINW)
        NCYC0 = ITER
        IF (IOUT.NE.3) GOTO 110
        CFAC = RES/(RESOLD+1.0D-40)
        RESOLD = RES
        IF (1.EQ.M) GOTO 150
        CALL RESID(2,RESCG,A,U,F,IA,JA,IW,IMIN,IMAX,IMINW)
        WRITE (IUP,9010) ITER,RESCG,RES,CFAC
        GOTO 110
150     WRITE (IUP,9020) ITER,RES,CFAC
110     IF (ICONV.EQ.1) GOTO 120
        EPSIL = EPSI
        IF (ICONV.NE.4) GOTO 115
        UMAX = 0.D0
        DO 160 I=IMIN(1),IMAX(1)
          UMAX = MAX(UMAX,ABS(U(I)))
160     CONTINUE
        EPSIL = EPSI*UMAX
115     IF (RES.LT.EPSIL) GOTO 170
120   CONTINUE
170   IF (IOUT.NE.3.AND.IOUT.NE.0) CALL RESID(1,RES,A,U,F,IA,JA,IW,
     +                                        IMIN,IMAX,IMINW)
      RETURN
C
9000  FORMAT (' CYCLE  0:',3X,'RES=',D9.3)
9005  FORMAT (/' ************* CYCLING..... *************'/)
9010  FORMAT ( ' CYCLE ',I2,':   RESCG=',D9.3,'   RES=',D9.3,
     +         '   CFAC=',D9.3)
9020  FORMAT ( ' CYCLE ',I2,':   RES=',D9.3,'   CFAC=',D9.3)
9030  FORMAT (' --------------------------------------------------',
     +        '----------')
9040  FORMAT (/' CYCLING BETWEEN GRIDS 1 AND',I3,':'/)
9050  FORMAT (' *** ERROR IN SOLVE: NDU TOO SMALL ***')
9060  FORMAT (' *** ERROR IN SOLVE: NDF TOO SMALL ***')
      END
C
C.......................................................................
C
C     CYC                                                    SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE CYC(L,NRDX,NRDTYP,NRDLEN,NRCX,NRUX,NRUTYP,NRULEN,
     +               IGAM,A,U,F,IA,JA,IW,IMIN,IMAX,IMINW,IMAXW,
     +               IFG,ICG,NSTCOL,NG,TIME,IROW0,M,IUM,IERR,ITER,
     +               NSC,NDA,NDJA,MDA,MDJA,MSEL,FAC,RESI,LEVELS)
C
C     PERFORMS ONE AMG CYCLE WITH GRID L AS FINEST GRID
C
C
      DOUBLE PRECISION A(*),U(*),F(*),RESI(*),RES,FAC
      REAL TIME(*),TNEW,TOLD
      INTEGER NRDTYP(*),NRUTYP(*),IA(*),JA(*),IW(*),IMIN(*),IMAX(*),
     +        IMINW(*),IMAXW(*),IFG(*),ICG(*),NSTCOL(*),NG(*)
      EXTERNAL COARSE,CTIME,INTA,PUTZ,RELX,RESID,VSCALE
      SAVE
C
C===> DURING FIRST CYCLE: INITIALIZE PARAMETERS CONTROLLING YALE-SMP
C     FACTORIZATION AND ADAPTIVE DETERMINATION OF COARSEST GRID:
C     IFAC=1: ON NEXT CALL OF YALE-SMP FACTORIZE MATRIX
C     IFI =1: ON FIRST RETURN TO NEXT TO COARSEST GRID AFTER COARSE
C             GRID SOLUTION COMPARE RESIDUAL WITH THE RESIDUAL ON THE
C             SAME GRID BEFORE COARSE GRID SOLUTION. IF REDUCTION OF
C             RESIDUAL NOT SATISFYING, REDUCE NUMBER OF GRIDS USED IN
C             CYCLING BY ONE AND REPETE THE PROCESS WITH THE NOW
C             COARSEST GRID.
C     NPTSF:  NUMBER OF GRID POINTS ON FINEST GRID USED IN CYCLE,
C             DIVIDED BY 10. ONLY GRIDS WITH LESS THEN NPTSF POINTS
C             ARE ALLOWED TO BECOME COARSEST GRID DURING ADAPTIVE
C             COARSE GRID DETERMINATION.
C
      IF (ITER.EQ.1) IFAC = 1
      IF (MSEL.NE.2) THEN
        IFI   = 0
        NPTSF = 0
      ELSE
        MINK  = 1000000
        IFI   = 1
        NPTSF = (IMAX(L)-IMIN(L)+1)/10
      ENDIF
      IF (L.LT.M) GOTO 100
C
C===> ONE GRID ONLY
C
      CALL COARSE(M,IFAC,NSC,NRCX,IUM,IERR,
     +            A,U,F,IA,JA,IW,IMIN,IMAX,IMINW,ICG,NSTCOL,TIME,
     +            NDA,NDJA,MDA,MDJA,IROW0)
      GOTO 1000
C
C===> MORE THAN ONE GRID
C
100   DO 110 K=L,M
        NG(K) = 0
110   CONTINUE
      IVSTAR = 3-IGAM
      K = L
C
C===> RELAX (DOWNWARDS)
C
150   DO 170 N=1,NRDX
        DO 160 NL=1,NRDLEN
          CALL RELX(K,NRDTYP(NL),A,U,F,IA,JA,IW,IMIN,IMAX,IMINW,ICG,
     +              NSTCOL,TIME)
160     CONTINUE
170   CONTINUE
      IF (IFI.NE.1.OR.IMAX(K)-IMIN(K).GE.NPTSF) GOTO 190
C
C===> MINK: LOWEST GRID NUMBER FOR WHICH RESIDUAL IS STORED DURING
C           FIRST DOWNWARDS RELAXATION
C
      IF (MINK.EQ.1000000) MINK = K
      CALL CTIME(TOLD)
      CALL RESID(K,RESI(K),A,U,F,IA,JA,IW,IMIN,IMAX,IMINW)
      CALL CTIME(TNEW)
      TIME(15) = TIME(15)+TNEW-TOLD
190   NG(K) = NG(K)+1
      K = K+1
195   CALL PUTZ(K,U,IMIN,IMAX,TIME)
      CALL RESC(K,A,U,F,IA,JA,IW,IMIN,IMAX,IMINW,IMAXW,IFG,TIME)
      IF (K.LT.M) GOTO 150
C
C===> SOLVE ON COARSEST GRID
C
      CALL COARSE(M,IFAC,NSC,NRCX,IUM,IERR,
     +            A,U,F,IA,JA,IW,IMIN,IMAX,IMINW,ICG,NSTCOL,TIME,
     +            NDA,NDJA,MDA,MDJA,IROW0)
C
C===> RELAX (UPWARDS)
C
200   CALL VSCALE(K,IVSTAR,A,U,F,IA,JA,IW,IMIN,IMAX,IMINW,TIME)
      K = K-1
      CALL INTA(K,A,U,IA,JA,IW,IMIN,IMAX,IMINW,IMAXW,IFG,TIME)
      DO 215 N=1,NRUX
        DO 210 NL=1,NRULEN
          CALL RELX(K,NRUTYP(NL),A,U,F,IA,JA,IW,IMIN,IMAX,IMINW,ICG,
     +              NSTCOL,TIME)
210     CONTINUE
215   CONTINUE
      IF (IFI.NE.1.OR.K.LT.MINK) GOTO 219
C
C===> ON FIRST RETURN TO NEXT TO COARSEST GRID COMPARE RESIDUAL WITH
C     THE PREVIOUS ONE
C
      CALL CTIME(TOLD)
      CALL RESID(K,RES,A,U,F,IA,JA,IW,IMIN,IMAX,IMINW)
      CALL CTIME(TNEW)
      TIME(15) = TIME(15)+TNEW-TOLD
C
C===> IF RESIDUAL REDUCTION SATISFYING: COARSE GRID ADAPTION FINISHED
C
      IF (RES.LT.RESI(K)*FAC) THEN
        IFI  = 0
      ELSE
        IF (NSC.EQ.2) LEVELS = K
        M    = K
        IFAC = 1
        GOTO 195
      ENDIF
219   IF (K.EQ.L) GOTO 1000
C
C===> GRID SWITCHING CORRESPONDING TO IGAM
C
      IF (IGAM.GE.3) GOTO 220
      IF (IGAM.GE.0) GOTO 200
      IF (K.EQ.L+1 .AND. NG(K).LT.ABS(IGAM)) GOTO 150
      GOTO 200
220   IF (NG(K).LT.2) GOTO 150
      IF (IGAM.EQ.4)  NG(K) = 0
      GOTO 200
C
C===> RETURN
C
1000  CALL NRMU(L,U,IMIN,IMAX,IROW0)
      END
C
C.......................................................................
C
C     COARSE                                                SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE COARSE(M,IFAC,NSC,NRCX,IUM,IERR,A,U,F,IA,JA,
     +                  IW,IMIN,IMAX,IMINW,ICG,NSTCOL,TIME,NDA,NDJA,
     +                  MDA,MDJA,IROW0)
C
C     SOLVES ON COARSEST GRID, EITHER WITH GAUSS-SEIDEL RELAXATION
C     (NSC=1) OR WITH THE YALE-SMP DIRECT SOLVER NDRV (NSC=2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION A(*),U(*),F(*)
      REAL TIME(*),TOLD,TNEW
      INTEGER IA(*),JA(*),IW(*),IMIN(*),IMAX(*),IMINW(*),ICG(*),
     +        NSTCOL(*),ESP,PATH,FLAG
      EXTERNAL CTIME,NDRV,RELX,RESID
C
C===> CONV: IF COARSE GRID SOLUTION IS DONE WITH GS-RELAXATION AND
C     NRCX=0, AS MANY GS-SWEEPS ARE PERFORMED AS ARE NECESSARY TO RE-
C     DUCE THE RESIDUAL BY THE FACTOR CONV
C
      PARAMETER (CONV=1.D-2)
      SAVE
C
      IF (NSC.NE.2) GOTO 400
C
C===> SOLUTION WITH YALE-SMP
C
      CALL CTIME(TOLD)
      ILO = IMIN(M)
      JLO = IA(ILO)
      IF (IFAC.EQ.1) THEN
C
C===> FIRST CALL ON GRID M, FIRST FACTORIZE MATRIX
C
        IHI = IMAX(M)
        JHI = IW(IMINW(M))-1
        NP = IHI-ILO+1
        IS = ILO-1
        JS = JLO-1
C
C===>   TEST OF AVAILABLE WORK SPACE
C
        IF (JHI+3*NP.GT.NDJA) THEN
          NSC = 1
          WRITE (IUM,9020)
          IERR = -3
          GOTO 400
        ENDIF
C
C===>   INITIALISATION OF YALE-SMP POINTER VECTORS
C
        DO 310 I=1,NP
          JA(JHI+I) = I
          JA(JHI+NP+I) = I
          JA(JHI+2*NP+I)= I
310     CONTINUE
        IF (IROW0.NE.1) THEN
C
C===>   COARSE GRID OPERATOR REGULAR, SHIFT CONTENTS OF POINTER
C       VECTORS IA AND JA
C
          DO 320 I=ILO,IHI
            IA(I) = IA(I)-JS
320       CONTINUE
          IAUX = IA(IHI+1)
          IA(IHI+1) = IW(IMINW(M))-JS
          DO 300 I=JLO,JHI
            JA(I) = JA(I)-IS
300       CONTINUE
          NPOINT = NP
        ELSE
C
C===>   COARSE GRID OPERATOR HAS ROWSUM ZERO, ELIMINATE LAST SOLUTION
C       COMPONENT BY SETTING U(IHI) TO ZERO AND CANCELLING THE COR-
C       RESPONDING ENTRIES IN A AND JA, RESPECTIVELY. THE CANCELLED
C       ENTRIES ARE STORED BEFORE THE LAST ROW IN A AND JA.
C
          IAUX = IA(IHI)
C
C===>     JPOS: POINTER TO POSITION IN A AND JA TO CONTAIN NEXT
C               ELIMINATED ENTRY
C
          JPOS = IAUX-1
          J = IA(ILO)
          DO 20 I=ILO,IHI-1
            IA(I) = J-JS
            JA(J) = I-IS
            J = J+1
13          IF (J.EQ.IA(I+1)) GOTO 20
            IF (JA(J).NE.IHI) THEN
              JA(J) = JA(J)-IS
              J = J+1
            ELSE
              AAUX = A(J)
              DO 11 JJ=J,JPOS-1
                A(JJ) = A(JJ+1)
                JA(JJ) = JA(JJ+1)
11            CONTINUE
              DO 12 II=I+1,IHI
                IA(II) = IA(II)-1
12            CONTINUE
              A(JPOS) = AAUX
              JA(JPOS) = I
              JPOS = JPOS-1
            ENDIF
            GOTO 13
20        CONTINUE
          IA(IHI) = IA(IHI)-JS
C
C===>     DECREASE NUMBER OF POINTS BY ONE AND SET LAST SOLUTION
C         COMPONENT TO ZERO
C
          NPOINT = NP-1
          U(IHI) = 0.D0
        ENDIF
        NSP = NDA-JHI
        PATH = 1
        CALL NDRV(NPOINT,JA(JHI+1),JA(JHI+NP+1),JA(JHI+2*NP+1),IA(ILO),
     +            JA(JLO),A(JLO),F(ILO),U(ILO),NSP,A(JHI+1),A(JHI+1),
     +            ESP,PATH,FLAG)
C
C===>   RESTORE PREVIOUS VALUES FOR IA, JA AND A, IN PARTICULAR PUT
C       BACK ELIMINATED MATRIX ENTRIES, IF ROWSUM ZERO
C
        DO 330 I=ILO,IHI
          IA(I) = IA(I)+JS
330     CONTINUE
        IF (IROW0.NE.1) THEN
          IA(IHI+1) = IAUX
          DO 335 I=JLO,JHI
            JA(I) = JA(I)+IS
335       CONTINUE
        ELSE
          DO 25 I=JLO,JPOS
            JA(I)=JA(I)+IS
25        CONTINUE
          DO 40 J=JPOS+1,IAUX-1
            AAUX = A(J)
            I = JA(J)
            DO 27 II=I+1,IHI
              IA(II) = IA(II)+1
27          CONTINUE
            DO 30 JJ=J,IA(I+1),-1
              A(JJ) = A(JJ-1)
              JA(JJ) = JA(JJ-1)
30          CONTINUE
            A(IA(I+1)-1) = AAUX
            JA(IA(I+1)-1) = IHI
40        CONTINUE
        ENDIF
C
C===>   IF AN ERROR OCCURED DURING EXECUTION OF NDRV, SOLVE WITH
C       GAUSS-SEIDEL RELAXATION
C
        IF (FLAG.NE.0) THEN
          NSC = 1
          IF (ESP.LT.0) THEN
            WRITE (IUM,9030)
            IERR = -1
          ELSE
            WRITE (IUM,9040) FLAG
            IERR = -32
          ENDIF
          GOTO 400
        ELSE
C
C===>     FACTORIZATION SUCCESSFULL, UPDATE LIMITS OF USED STORAGE
C
          MDA = MAX(MDA,NDA-ESP)
          MDJA = MAX(MDJA,JHI+3*NP)
          IFAC = 0
        ENDIF
      ELSE
C
C===>   FACTORIZATION ALLREADY DONE
C
        IF (IROW0.NE.1) THEN
          NPOINT = NP
        ELSE
          NPOINT = NP-1
          U(IHI) = 0.D0
        ENDIF
        PATH = 3
        CALL NDRV(NPOINT,JA(JHI+1),JA(JHI+NP+1),JA(JHI+2*NP+1),IA(ILO),
     +            JA(JLO),A(JLO),F(ILO),U(ILO),NSP,A(JHI+1),A(JHI+1),
     +            ESP,PATH,FLAG)
      ENDIF
C
C===> UPDATE TIME COUNTER
C
      CALL CTIME(TNEW)
      TIME(17) = TIME(17)+TNEW-TOLD
      GOTO 190
C
C===> SOLUTION WITH GAUSS-SEIDEL RELAXATION
C
400   IF (NRCX.NE.0) THEN
        DO 180 ITER=1,NRCX
          CALL RELX(M,2,A,U,F,IA,JA,IW,IMIN,IMAX,IMINW,ICG,NSTCOL,TIME)
180     CONTINUE
      ELSE
C
C===> IF NRCX=0: REDUCE RESIDUAL ON COARSEST GRID BY FACTOR CONV
C                (IF NOT YET IN THE RANGE OF THE TRUNCATION ERROR)
C
        CALL CTIME(TOLD)
C
C===>   CALCULATE SUPREMUM NORM OF RIGHT HAND SIDE
C
        FMAX = 0.D0
        DO 181 I=IMIN(M),IMAX(M)
          FMAX = MAX(FMAX,ABS(F(I)))
181     CONTINUE
        CALL RESID(M,RESOLD,A,U,F,IA,JA,IW,IMIN,IMAX,IMINW)
        CALL CTIME(TNEW)
        TIME(15) = TIME(15)+TNEW-TOLD
        RESOLD = MAX(RESOLD*CONV,FMAX*1.D-12)
        DO 187 I=1,10
          DO 185 J=1,10
            CALL RELX(M,2,A,U,F,IA,JA,IW,IMIN,IMAX,IMINW,ICG,NSTCOL,
     +                TIME)
185       CONTINUE
          CALL CTIME(TOLD)
          CALL RESID(M,RESNEW,A,U,F,IA,JA,IW,IMIN,IMAX,IMINW)
          CALL CTIME(TNEW)
          TIME(15) = TIME(15)+TNEW-TOLD
          IF (RESNEW.LE.RESOLD) GOTO 190
187     CONTINUE
      ENDIF
C
C===> COMPUTE RESIDUAL AFTER SOLUTION ON GRID M
C
190   RETURN
C
9020  FORMAT (' --- WARNG IN COARSE: NO YALE-SMP BECAUSE NDJA TOO ',
     +        'SMALL')
9030  FORMAT (' --- WARNG IN COARSE: NO YALE-SMP BECAUSE NDA TOO SMALL')
9040  FORMAT (' --- WARNG IN COARSE: NO YALE-SMP BECAUSE OF ERROR IN ',
     +        'FACTORIZATION,'/'     CODE=',I8)
      END
C
C.......................................................................
C
C     RELX                                                  SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE RELX(K,IREL,A,U,F,IA,JA,IW,IMIN,IMAX,IMINW,ICG,NSTCOL,
     +                TIME)
C
C     PERFORMS ONE (PARTIAL) GAUSS-SEIDEL SWEEP ON GRIK K:
C
C     IREL = 1 :   PARTIAL GAUSS-SEIDEL SWEEP (ONLY F-POINTS)
C          = 2 :   FULL GAUSS-SEIDEL SWEEP (ALL POINTS)
C          = 3 :   PARTIAL GAUSS-SEIDEL SWEEP (ONLY C-POINTS)
C          = 4 :   FULL SWEEP: FF -- C -- COLORS (HIGHEST FIRST)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION A(*),U(*),F(*)
      REAL TOLD,TNEW,TIME(*)
      INTEGER IA(*),JA(*),IW(*),IMIN(*),IMAX(*),IMINW(*),ICG(*),
     +        NSTCOL(*)
      EXTERNAL CTIME
      SAVE
C
      CALL CTIME(TOLD)
      IAUX = IA(IMAX(K)+1)
      IA(IMAX(K)+1) = IW(IMINW(K))
      GOTO (100,200,300,400), IREL
      GOTO 200
C
C===> F-RELAXATION
C
100   DO 120 I=IMIN(K),IMAX(K)
        IF (ICG(I).GT.0) GOTO 120
        S = F(I)
        DO 110 J=IA(I)+1,IA(I+1)-1
          S = S-A(J)*U(JA(J))
110     CONTINUE
        U(I) = S/A(IA(I))
120   CONTINUE
      GOTO 1000
C
C===> FULL GS RELAXATION
C
200   DO 220 I=IMIN(K),IMAX(K)
        S = F(I)
        DO 210 J=IA(I)+1,IA(I+1)-1
          S = S-A(J)*U(JA(J))
210     CONTINUE
        U(I) = S/A(IA(I))
220   CONTINUE
      GOTO 1000
C
C===> C-RELAXATION
C
300   DO 320 I=IMIN(K),IMAX(K)
        IF (ICG(I).LE.0) GOTO 320
        S = F(I)
        DO 310 J=IA(I)+1,IA(I+1)-1
          S = S-A(J)*U(JA(J))
310     CONTINUE
        U(I) = S/A(IA(I))
320   CONTINUE
      GOTO 1000
C
C===> FF-RELAXATION
C
400   DO 420 I=IMIN(K),IMAX(K)
        IF (ICG(I).NE.0) GOTO 420
        S = F(I)
        DO 410 J=IA(I)+1,IA(I+1)-1
          S = S-A(J)*U(JA(J))
410     CONTINUE
        U(I) = S/A(IA(I))
420   CONTINUE
C
C===> C-RELAXATION
C
      DO 440 I=IMIN(K),IMAX(K)
        IF (ICG(I).LE.0) GOTO 440
        S = F(I)
        DO 430 J=IA(I)+1,IA(I+1)-1
          S = S-A(J)*U(JA(J))
430     CONTINUE
        U(I) = S/A(IA(I))
440   CONTINUE
C
C===> MORE-COLOR RELAXATION
C
      I = NSTCOL(K)
470   IF (I.GE.100000000) GOTO 1000
      S = F(I)
      DO 480 J=IA(I)+1,IA(I+1)-1
        S = S-A(J)*U(JA(J))
480   CONTINUE
      U(I) = S/A(IA(I))
      I = -ICG(I)
      GOTO 470
C
1000  CALL CTIME(TNEW)
      TIME(13) = TIME(13)+TNEW-TOLD
      IA(IMAX(K)+1) = IAUX
      RETURN
      END
C
C.......................................................................
C
C     VSCALE                                                SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE VSCALE(K,IVSTAR,A,U,F,IA,JA,IW,IMIN,IMAX,IMINW,TIME)
C
C     SCALES ACTUAL APPROXIMATE SOLUTION ON GRID K (V*-CYCLE); SCALING
C     IS DONE SUCH THAT ENERGY NORM BECOMES MINIMAL
C
C     NOTE: THIS SCALING MAKES SENSE ONLY FOR SYMMETRIC PROBLEMS
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION A(*),U(*),F(*)
      REAL TOLD,TNEW,TIME(*)
      INTEGER IA(*),JA(*),IW(*),IMIN(*),IMAX(*),IMINW(*)
      EXTERNAL CTIME
      SAVE
C
C===> COMPUTATION OF SCALING FACTOR "FAC"
C
      IF (IVSTAR.NE.1) RETURN
      CALL CTIME(TOLD)
      IAUX = IA(IMAX(K)+1)
      IA(IMAX(K)+1) = IW(IMINW(K))
      S1 = 0.0D0
      S2 = 0.0D0
      DO 10 I=IMIN(K),IMAX(K)
        SA = 0.0D0
        DO 20 J=IA(I),IA(I+1)-1
          SA = SA+A(J)*U(JA(J))
20      CONTINUE
        S1 = S1+U(I)*F(I)
        S2 = S2+U(I)*SA
10    CONTINUE
      FAC = 1.0D0
      IF (S2.NE.0.0D0) FAC = S1/S2
C
C===> SCALING
C
      DO 30 I=IMIN(K),IMAX(K)
        U(I) = U(I)*FAC
30    CONTINUE
      IA(IMAX(K)+1) = IAUX
      CALL CTIME(TNEW)
      TIME(14) = TIME(14)+TNEW-TOLD
      RETURN
      END
C
C.......................................................................
C
C     INTA                                                  SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE INTA(KF,A,U,IA,JA,IW,IMIN,IMAX,IMINW,IMAXW,IFG,TIME)
C
C     INTERPOLATES CORRECTION FROM GRID KF+1 TO GRID KF
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION A(*),U(*)
      REAL TOLD,TNEW,TIME(*)
      INTEGER IA(*),JA(*),IW(*),IMIN(*),IMAX(*),IMINW(*),IMAXW(*),IFG(*)
      EXTERNAL CTIME
      SAVE
C
C===> C->C CONTRIBUTIONS
C
      CALL CTIME(TOLD)
      DO 50 IC=IMIN(KF+1),IMAX(KF+1)
        IF = IFG(IC)
        U(IF) = U(IF)+U(IC)
50    CONTINUE
C
C===> C->F CONTRIBUTIONS
C
      IAUX = IW(IMAXW(KF)+1)
      IW(IMAXW(KF)+1) = IA(IMIN(KF+1))
      DO 100 I=IMINW(KF),IMAXW(KF)
        IF = IFG(I)
        DO 150 J=IW(I),IW(I+1)-1
          U(IF) = U(IF)+A(J)*U(JA(J))
150     CONTINUE
100   CONTINUE
      IW(IMAXW(KF)+1) = IAUX
      CALL CTIME(TNEW)
      TIME(11) = TIME(11)+TNEW-TOLD
      RETURN
      END
C
C.......................................................................
C
C     RESC                                                  SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE RESC(KC,A,U,F,IA,JA,IW,IMIN,IMAX,IMINW,IMAXW,IFG,TIME)
C
C     RESTRICTS RESIDUALS FROM GRID KC-1 TO GRID KC
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION A(*),U(*),F(*)
      REAL TOLD,TNEW,TIME(*)
      INTEGER IA(*),JA(*),IW(*),IMIN(*),IMAX(*),IMINW(*),IMAXW(*),IFG(*)
      EXTERNAL CTIME
      SAVE
C
C===> TRANSFER OF C-POINT DEFECTS
C
      CALL CTIME(TOLD)
      IAUX = IA(IMAX(KC-1)+1)
      IA(IMAX(KC-1)+1) = IW(IMINW(KC-1))
      IAUX1 = IW(IMAXW(KC-1)+1)
      IW(IMAXW(KC-1)+1) = IAUX
      DO 100 IC=IMIN(KC),IMAX(KC)
        IF = IFG(IC)
        D = F(IF)
        DO 80 J=IA(IF),IA(IF+1)-1
          D = D-A(J)*U(JA(J))
80      CONTINUE
        F(IC) = D
100   CONTINUE
C
C===> TRANSFER OF F-POINT DEFECTS
C
      DO 200 I=IMINW(KC-1),IMAXW(KC-1)
        IF = IFG(I)
        D = F(IF)
        DO 20 J=IA(IF),IA(IF+1)-1
          D = D-A(J)*U(JA(J))
20      CONTINUE
        DO 250 J=IW(I),IW(I+1)-1
          F(JA(J)) = F(JA(J))+A(J)*D
250     CONTINUE
200   CONTINUE
      IA(IMAX(KC-1)+1) = IAUX
      IW(IMAXW(KC-1)+1) = IAUX1
      CALL CTIME(TNEW)
      TIME(12) = TIME(12)+TNEW-TOLD
      RETURN
      END
C
C.......................................................................
C
C     PUTZ                                                   SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE PUTZ(K,U,IMIN,IMAX,TIME)
C
C     PUTS ZERO TO U-VALUES OF GRID K
C
      DOUBLE PRECISION U(*)
      REAL TOLD,TNEW,TIME(*)
      INTEGER IMIN(*),IMAX(*)
      EXTERNAL CTIME
      SAVE
C
      CALL CTIME(TOLD)
      DO 10 I=IMIN(K),IMAX(K)
        U(I) = 0.0D0
10    CONTINUE
      CALL CTIME(TNEW)
      TIME(15) = TIME(15)+TNEW-TOLD
      RETURN
      END
C
C.......................................................................
C
C     FIRST                                                 SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE FIRST(IFIRST,U,IMIN,IMAX,IARR,IROW0)
C
C     PUTS A FIRST APPROXIMATION TO FINEST GRID
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION U(*)
      INTEGER IMIN(*),IMAX(*),IARR(*)
      EXTERNAL IDEC,RANDOM
      SAVE
C
      IF (IFIRST.NE.0) THEN
        CALL IDEC(IFIRST,3,NDIGIT,IARR)
        IFRST = IARR(2)
      ELSE
        IFRST = 3
      ENDIF
C
      GOTO (100,200,300), IFRST
      RETURN
C
100   DO 110 I=IMIN(1),IMAX(1)
        U(I) = 0.0D0
110   CONTINUE
      RETURN
C
200   DO 210 I=IMIN(1),IMAX(1)
        U(I) = 1.0D0
210   CONTINUE
      IF (IROW0.LT.2) U(IMAX(1)) = 0.0D0
      RETURN
C
300   IF (IARR(3)*IFIRST.EQ.0) GOTO 350
      S = DBLE(IARR(3))
      DO 310 I=1,10
         S = S*0.1D0
         IF (S.LT.1.0D0) GOTO 370
310   CONTINUE
350   S = 0.72815D0
370   DO 390 I=IMIN(1),IMAX(1)
        U(I) = RANDOM(S)
390   CONTINUE
      IF (IROW0.LT.2) U(IMAX(1)) = 0.0D0
      RETURN
      END
C
C.......................................................................
C
C     INJF                                                   SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE INJF(KC,F,IMIN,IMAX,IFG)
C
C     INJECTS F-VALUES FROM GRID KC-1 TO GRID KC
C
      DOUBLE PRECISION F(*)
      INTEGER IMIN(*),IMAX(*),IFG(*)
      SAVE
      DO 10 IC=IMIN(KC),IMAX(KC)
        F(IC) = F(IFG(IC))
10    CONTINUE
      RETURN
      END
C
C.......................................................................
C
C     INJU                                                   SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE INJU(KC,U,IMIN,IMAX,IFG)
C
C     INJECTS U-VALUES FROM GRID KC-1 TO GRID KC
C
      DOUBLE PRECISION U(*)
      INTEGER IMIN(*),IMAX(*),IFG(*)
      SAVE
      DO 10 IC=IMIN(KC),IMAX(KC)
        U(IC) = U(IFG(IC))
10    CONTINUE
      RETURN
      END
C
C.......................................................................
C
C     NRMU                                                SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE NRMU(K,U,IMIN,IMAX,IROW0)
C
C     NORMALIZES U ON GRID K IF ROWSUM=0 (LAST COMPONENT =0)
C
      DOUBLE PRECISION U(*),FAC
      INTEGER IMIN(*),IMAX(*)
      SAVE
C
      IF (IROW0.GT.1) RETURN
      FAC = U(IMAX(K))
      DO 10 I=IMIN(K),IMAX(K)
        U(I) = U(I)-FAC
10    CONTINUE
      RETURN
      END
C
C.......................................................................
C
C     RESID                                                 SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE RESID(K,RES,A,U,F,IA,JA,IW,IMIN,IMAX,IMINW)
C
C     COMPUTES L2-NORM OF RESIDUAL ON GRID K
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION A(*),U(*),F(*)
      INTEGER IA(*),JA(*),IW(*),IMIN(*),IMAX(*),IMINW(*)
      SAVE
C
      RES = 0.0D0
      IAUX = IA(IMAX(K)+1)
      IA(IMAX(K)+1) = IW(IMINW(K))
      DO 20 I=IMIN(K),IMAX(K)
        S = F(I)
        DO 10 J=IA(I),IA(I+1)-1
          S = S-A(J)*U(JA(J))
10      CONTINUE
        RES = RES+S*S
20    CONTINUE
      IA(IMAX(K)+1) = IAUX
      RES = SQRT(RES)
      RETURN
      END
C
C.......................................................................
C
C     CG                                                 SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE CG(K,ICGR,ITER,A,U,F,IA,JA,IW,IMIN,IMAX,IMINW,M,TIME,
     +              IERR,IUM)
C
C     PERFORMS ONE STEP OF PRECONDITIONED CONJUGATE GRADIENT
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION A(*),U(*),F(*)
      REAL TOLD,TNEW,TIME(*)
      INTEGER IA(*),JA(*),IW(*),IMIN(*),IMAX(*),IMINW(*)
      EXTERNAL CGALF,CGEPS,CTIME
      SAVE
C
      IF (ICGR.EQ.0) RETURN
C
C===> COMPUTE MOST RECENT MG CORRECTION
C
      CALL CTIME(TOLD)
      NNU = IMAX(K)-IMIN(K)+1
      ISHIFT = IMAX(M)+1-IMIN(K)
      DO 10 I=IMIN(K),IMAX(K)
        U(I) = U(I)-U(I+ISHIFT)
10    CONTINUE
      IF (ICGR.EQ.2 .AND. ITER.GT.1) GOTO 100
C
C===> FIRST CG STEP
C
      DO 20 I=IMIN(K),IMAX(K)
        F(I+ISHIFT) = U(I)
20    CONTINUE
      GOTO 200
C
C===> NEXT CG STEPS (IF ICGR=2 ONLY)
C
100   ALF = CGALF(K,S2,A,U,F,IA,JA,IW,IMIN,IMAX,IMINW,M)
      DO 110 I=IMIN(K),IMAX(K)
        F(I+ISHIFT) = U(I)+ALF*F(I+ISHIFT)
110   CONTINUE
200   EPS = CGEPS(K,S2,A,U,F,IA,JA,IW,IMIN,IMAX,IMINW,M,IERR,IUM)
      IF (IERR.GT.0) RETURN
      DO 210 I=IMIN(K),IMAX(K)
        U(I) = U(I+ISHIFT)+EPS*F(I+ISHIFT)
210   CONTINUE
      CALL CTIME(TNEW)
      TIME(16) = TIME(16)+TNEW-TOLD
      RETURN
      END
C
C.......................................................................
C
C     USAVE                                            SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE USAVE(K,ICGR,U,IMIN,IMAX,NDU,M,TIME,IERR,IUM,MDU,NDF,
     +                 MDF)
C
C     MAKES A BACK-UP OF THE CURRENT APPROXIMATION ON LEVEL K
C     IF ICGR.NE.0.
C
      DOUBLE PRECISION U(*)
      REAL TIME(*)
      INTEGER IMIN(*),IMAX(*)
      EXTERNAL CTIME
      SAVE
C
      IF (ICGR.EQ.0) RETURN
      ISHIFT = IMAX(M)+1-IMIN(K)
      MDU = MAX(MDU,IMAX(K)+ISHIFT)
      MDF = MAX(MDF,IMAX(K)+ISHIFT)
      IF (MDU.LE.NDU.AND.MDF.LE.NDF) GOTO 10
      WRITE (IUM,9000) MDU,MDF
      IERR = -4
      ICGR = 0
      RETURN
C
10    CALL CTIME(TOLD)
      DO 20 I=IMIN(K),IMAX(K)
        U(I+ISHIFT) = U(I)
20    CONTINUE
      CALL CTIME(TNEW)
      TIME(16) = TIME(16)+TNEW-TOLD
      RETURN
C
9000  FORMAT (' --- WARNG IN USAVE: NO CG BECAUSE OF STORAGE ---'/
     +        '     REQUIRED: NDU =',    I9                      /
     +        '               NDF =',    I9                      )
      END
C
C.......................................................................
C
C     CGEPS                                            FUNCTION
C
C.......................................................................
C
      DOUBLE PRECISION FUNCTION CGEPS(K,S2,A,U,F,IA,JA,IW,
     +                                IMIN,IMAX,IMINW,M,IERR,IUM)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION A(*),U(*),F(*)
      INTEGER IMIN(*),IMAX(*),IMINW(*),IA(*),IW(*),JA(*)
      SAVE
C
      ISHIFT = IMAX(M)+1-IMIN(K)
      S1 = 0.0D0
      S2 = 0.0D0
      IAUX = IA(IMAX(K)+1)
      IA(IMAX(K)+1) = IW(IMINW(K))
      DO 50 I=IMIN(K),IMAX(K)
        SR = F(I)
        SP = 0.0D0
        DO 40 J=IA(I),IA(I+1)-1
          SR = SR-A(J)*U(JA(J)+ISHIFT)
          SP = SP+A(J)*F(JA(J)+ISHIFT)
40      CONTINUE
        S1 = S1+SR*F(I+ISHIFT)
        S2 = S2+SP*F(I+ISHIFT)
50    CONTINUE
      IA(IMAX(K)+1) = IAUX
      IF (S2.EQ.0.0D0) GOTO 100
      CGEPS = S1/S2
      RETURN
C
C===> ERROR EXIT
C
100   WRITE (IUM,9000)
      IERR = 31
      RETURN
9000  FORMAT (' *** ERROR IN CGEPS: CG CORRECTION NOT DEFINED ***')
      END
C
C.......................................................................
C
C     CGALF                                            FUNCTION
C
C.......................................................................
C
      DOUBLE PRECISION FUNCTION CGALF(K,S2,A,U,F,IA,JA,IW,
     +                                IMIN,IMAX,IMINW,M)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION A(*),U(*),F(*)
      INTEGER IA(*),JA(*),IW(*),IMIN(*),IMAX(*),IMINW(*)
      SAVE
C
      ISHIFT=IMAX(M)+1-IMIN(K)
      S1 = 0.0D0
      IAUX = IA(IMAX(K)+1)
      IA(IMAX(K)+1) = IW(IMINW(K))
      DO 50 I=IMIN(K),IMAX(K)
        SR = 0.0D0
        DO 40 J=IA(I),IA(I+1)-1
          SR = SR+A(J)*U(JA(J))
40      CONTINUE
        S1 = S1+SR*F(I+ISHIFT)
50    CONTINUE
      CGALF = -S1/S2
      IA(IMAX(K)+1) = IAUX
      RETURN
      END
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     OUTPUT OF STATISTICAL INFORMATION ON CP-TIMES AND DIMENSIONING
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C.......................................................................
C
C     WRKCNT                                           SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE WRKCNT(IOUT,IA,IW,IMIN,IMAX,IMINW,LEVELS,TIME,NCYC0,
     +                  IUP,MDA,MDIA,MDJA,MDU,MDF,MDIG,RES0,RES)
C
C     RESIDUALS / CP-TIMES / COMPLEXITY / DIMENSIONING
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      REAL TIME(*),T(10),SUM1,SUM2
      INTEGER IA(*),IW(*),IMIN(*),IMAX(*),IMINW(*)
      SAVE
C
      IF (IOUT.LT.1) RETURN
C
C===> RESIDUALS / CONVERGENCE
C
      IF (NCYC0.GT.0) THEN
        CFAC = RES/(RES0+1.D-40)
        WRITE (IUP,9110) RES0,RES,CFAC
        CFPC = CFAC**(1.D0/DBLE(NCYC0))
        WRITE (IUP,9120) CFPC
      ENDIF
      IF (IOUT.LE.1) RETURN
      WRITE (IUP,9000)
      NNU = IMAX(1)
C
C===> COMPUTING TIMES
C
      SUM1 = 0.0
      SUM2 = 0.0
      DO 10 I=1,10
        T(I) = 0.0
        IF (NCYC0.GT.0) T(I) = TIME(I+10)/REAL(NCYC0)
        SUM1 = SUM1+TIME(I)
        SUM2 = SUM2+T(I)
10    CONTINUE
C
      WRITE (IUP,9020)
      WRITE (IUP,9030)
      WRITE (IUP,9040) (TIME(I),T(I),I=1,8)
      WRITE (IUP,9030)
      WRITE (IUP,9050) SUM1,SUM2
      WRITE (IUP,9030)
C
C===> SPACE OCCUPIED BY OPERATORS A(1) - A(LEVELS)
C
      IDIMA = 0
      DO 20 K=1,LEVELS
        IDIMA = IDIMA+IW(IMINW(K))-IA(IMIN(K))
20    CONTINUE
C
C===> THEORETICAL MINIMAL SPACE REQUIREMENTS
C
      IF (LEVELS.LT.2) RETURN
      MDTA  = IW(IMINW(LEVELS))-1
      MDTJA = MDTA
      MDTIA = IMAX(LEVELS)+1
      MDTU  = IMAX(LEVELS)
      MDTF  = MDTU
      MDTIG = 2*MDTU+NNU
      WRITE (IUP,9100) MDA,MDTA,MDJA,MDTJA,MDIA,MDTIA,MDU,MDTU,
     +                 MDF,MDTF,MDIG,MDTIG
C
C===> COMPLEXITIES
C
      SCMPLX = DBLE(2*(MDA+MDU+MDF)+MDJA+MDIA+MDIG)/
     +         DBLE(1+5*NNU+3*(IW(IMINW(1))-IA(IMIN(1))))
      TCMPLX = DBLE(2*(MDTA+MDTU+MDTF)+MDTJA+MDTIA+MDTIG)/
     +         DBLE(1+5*NNU+3*(IW(IMINW(1))-IA(IMIN(1))))
C
      ACMPLX = DBLE(IDIMA)/DBLE(IW(1)-1)
      OCMPLX = DBLE(MDTU )/DBLE(NNU)
      WRITE (IUP,9080) ACMPLX,OCMPLX,SCMPLX,TCMPLX
      RETURN
C
9000  FORMAT (//' ************** WORK COUNT ***************'/)
9020  FORMAT (  '   PREP       SEC       SOL      SEC/CYCLE')
9030  FORMAT (  ' -----------------------------------------')
9040  FORMAT (  ' 1 RWSRT   ',F7.2,'   11 INTADD   ',F7.2/
     +          ' 2 PRE-COL ',F7.2,'   12 RESCAL   ',F7.2/
     +          ' 3 CHK-COL ',F7.2,'   13 RELAX    ',F7.2/
     +          ' 4 INTERPOL',F7.2,'   14 V-*      ',F7.2/
     +          ' 5 RESTRICT',F7.2,'   15 OTHERS   ',F7.2/
     +          ' 6 OPDFN   ',F7.2,'   16 CONJ-GRAD',F7.2/
     +          ' 7 TRUNC   ',F7.2,'   17 YALE-SMP ',F7.2/
     +          ' 8 OTHERS  ',F7.2,'   18 ------   ',F7.2)
9050  FORMAT (  '   SUM     ',F7.2,'      SUM      ',F7.2)
9080  FORMAT (/' ******************* COMPLEXITIES ********************'/
     +         ' SPACE OCCUPIED BY ALL OPERATORS / SPACE OF OPERATOR  '/
     +         ' ON THE FINEST GRID   = ',F8.2,'   (A-COMPLEXITY)     '/
     +         ' TOTAL NUMBER OF GRID POINTS / NUMBER OF POINTS IN    '/
     +         ' THE  FINEST  GRID    = ',F8.2,'   (O-COMPLEXITY)     '/
     +         ' TOTAL SPACE USED BY AMG1R5 / SPACE OCCUPIED BY USER- '/
     +         ' DEFINED  PROBLEM     = ',F8.2,'   (S-COMPLEXITY)     '/
     +         ' SPACE USED DURING SOLUTION PHASE / SPACE OCCUPIED BY '/
     +         ' USER-DEFINED PROBLEM = ',F8.2                         /
     +         ' *****************************************************')
9100  FORMAT (//' ********* SPACE REQUIREMENTS *********'//
     +          ' VECTOR      NEEDED      THEOR. MINIMUM'/
     +          ' --------------------------------------'/
     +          '   A ',        I14   ,2X,     I16       /
     +          '   JA',        I14   ,2X,     I16       /
     +          '   IA',        I14   ,2X,     I16       /
     +          '   U ',        I14   ,2X,     I16       /
     +          '   F ',        I14   ,2X,     I16       /
     +          '   IG',        I14   ,2X,     I16       /
     +          ' --------------------------------------'/)
9110  FORMAT (//' **************** CONVERGENCE *****************'/
     +          ' L2-NORM OF RESIDUAL BEFORE CYCLING =',D10.3/
     +          ' L2-NORM OF RESIDUAL AFTER  CYCLING =',D10.3/
     +          ' CONVERGENCE FACTOR                 =',D10.3)
9120  FORMAT (  ' CONVERGENCE FACTOR PER CYCLE       =',D10.3)
      END
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     AMG1R5 AUXILIARY PROGRAMS
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C.......................................................................
C
C     IDEC                                                  SUBROUTINE
C
C.......................................................................
C
      SUBROUTINE IDEC(INTO,NNUM,NDIGIT,IARR)
C
C     DECOMPOSE NON-NEGATIVE INTEGER INTO INTO NNUM INTEGERS
C
C     INPUT:  INTO   - INTEGER (0.LE. INTO .LE.999999999)
C             NNUM   - INTEGER (1.LE. NNUM .LE.9); NUMBER OF INTEGERS
C                      TO BE RETURNED ON ARRAY IARR (SEE BELOW)
C
C     OUTPUT: NDIGIT - INTEGER; NUMBER OF DIGITS OF INTO
C             IARR   - INTEGER-ARRAY OF LENGTH 10:
C                      IARR(1)        = FIRST      DIGIT OF INTO,
C                      IARR(2)        = SECOND     DIGIT OF INTO, .....
C                      IARR(NNUM-1)   = (NNUM-1)ST DIGIT OF INTO,
C                      IARR(NNUM)     = REST OF INTO
C                      IF NNUM > NDIGIT, THE CORRESPONDING COMPONENTS
C                      OF IARR ARE PUT TO ZERO.
C
C     WARNING: BE SURE THAT YOUT COMPUTER CAN STORE NNUM DIGITS ON AN
C              INTEGER VARIABLE.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION IARR(*)
      DATA EPS /0.5D0/
      SAVE
C
      NREST = INTO
      NDIGIT = 1+INT(LOG10(EPS+DBLE(INTO)))
      IF (NNUM.GE.NDIGIT) GOTO 20
      IQ = 10**(NDIGIT-NNUM+1)
      IARR(NNUM) = NREST-(NREST/IQ)*IQ
      NREST = NREST/IQ
      DO 10 I=NNUM-1,1,-1
        IARR(I) = NREST-(NREST/10)*10
        NREST = NREST/10
10    CONTINUE
      RETURN
C
20    DO 30 I=NDIGIT,1,-1
        IARR(I) = NREST-(NREST/10)*10
        NREST = NREST/10
30    CONTINUE
      DO 40 I=NDIGIT+1,NNUM
        IARR(I) = 0
40    CONTINUE
      RETURN
      END
C
C.......................................................................
C
C     RANDOM                                                FUNCTION
C
C.......................................................................
C
      DOUBLE PRECISION FUNCTION RANDOM(S)
C
C     FUNCTION TO CREATE "RANDOM" SEQUENCE OF NUMBERS BETWEEN 0 AND 0.1
C
C     INPUT:   S      - NUMBER BETWEEN 0 AND 0.1
C
C     OUTPUT:  RANDOM - NUMBER BETWEEN 0 AND 0.1
C              S      - S=RANDOM
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SAVE
      RANDOM=100.0D0*DEXP(S)
      RANDOM=RANDOM-DBLE(INT(RANDOM))
      S=RANDOM
      RETURN
      END
c***********************************************************************
C $LARGE
C $NOFLOATCALLS
C                            APPENDIX 1                          1/15/81
C
C        SUBROUTINES FOR SOLVING SPARSE NONSYMMETRIC SYSTEMS
C        OF LINEAR EQUATIONS  (UNCOMPRESSED POINTER STORAGE)
C
C        REAL*8 VERSION. NOTE: THE ORIGINAL SUBROUTINES
C
C            NDRV, NSF, NNF, NNS AND NNT
C
C        HAVE BEEN RENAMED TO
C
C            YALE8, NSF8, NNF8, NNS8 AND NNT8, RESPECTIVELY.
C
C*** SUBROUTINE YALE8 (OLD NAME: NDRV)
C*** DRIVER FOR SUBROUTINES FOR SOLVING SPARSE NONSYMMETRIC SYSTEMS OF
C       LINEAR EQUATIONS (UNCOMPRESSED POINTER STORAGE)
C
C       SUBROUTINE  NDRV  (= OLD NAME)
C       SUBROUTINE  YALE8 (= NEW NAME)
        SUBROUTINE  NDRV
     *     (N, R,C,IC, IA,JA,A, B, Z, NSP,ISP,RSP,ESP, PATH, FLAG)
C
C    PARAMETERS
C    CLASS ABBREVIATIONS ARE --
C       N - INTEGER VARIABLE
C       F - REAL VARIABLE
C       V - SUPPLIES A VALUE TO THE DRIVER
C       R - RETURNS A RESULT FROM THE DRIVER
C       I - USED INTERNALLY BY THE DRIVER
C       A - ARRAY
C
C CLASS   PARAMETER
C ------+----------
C
C         THE NONZERO ENTRIES OF THE COEFFICIENT MATRIX M ARE STORED
C    ROW-BY-ROW IN THE ARRAY A.  TO IDENTIFY THE INDIVIDUAL NONZERO
C    ENTRIES IN EACH ROW, WE NEED TO KNOW IN WHICH COLUMN EACH ENTRY
C    LIES.  THE COLUMN INDICES WHICH CORRESPOND TO THE NONZERO ENTRIES
C    OF M ARE STORED IN THE ARRAY JA;  I.E., IF  A(K) = M(I,J),  THEN
C    JA(K) = J.  IN ADDITION, WE NEED TO KNOW WHERE EACH ROW STARTS AND
C    HOW LONG IT IS.  THE INDEX POSITIONS IN JA AND A WHERE THE ROWS OF
C    M BEGIN ARE STORED IN THE ARRAY IA;  I.E., IF M(I,J) IS THE FIRST
C    NONZERO ENTRY (STORED) IN THE I-TH ROW AND A(K) = M(I,J),  THEN
C    IA(I) = K.  MOREOVER, THE INDEX IN JA AND A OF THE FIRST LOCATION
C    FOLLOWING THE LAST ELEMENT IN THE LAST ROW IS STORED IN IA(N+1).
C    THUS, THE NUMBER OF ENTRIES IN THE I-TH ROW IS GIVEN BY
C    IA(I+1) - IA(I),  THE NONZERO ENTRIES OF THE I-TH ROW ARE STORED
C    CONSECUTIVELY IN
C            A(IA(I)),  A(IA(I)+1),  ..., A(IA(I+1)-1),
C    AND THE CORRESPONDING COLUMN INDICES ARE STORED CONSECUTIVELY IN
C            JA(IA(I)), JA(IA(I)+1), ..., JA(IA(I+1)-1).
C    FOR EXAMPLE, THE 5 BY 5 MATRIX
C                ( 1. 0. 2. 0. 0.)
C                ( 0. 3. 0. 0. 0.)
C            M = ( 0. 4. 5. 6. 0.)
C                ( 0. 0. 0. 7. 0.)
C                ( 0. 0. 0. 8. 9.)
C    WOULD BE STORED AS
C                 1  2  3  4  5  6  7  8  9
C            ---+--------------------------
C            IA   1  3  4  7  8 10
C            JA   1  3  2  2  3  4  4  4  5
C             A   1. 2. 3. 4. 5. 6. 7. 8. 9.         .
C
C NV      N     - NUMBER OF VARIABLES/EQUATIONS.
C FVA     A     - NONZERO ENTRIES OF THE COEFFICIENT MATRIX M, STORED
C                   BY ROWS.
C                   SIZE = NUMBER OF NONZERO ENTRIES IN M.
C NVA     IA    - POINTERS TO DELIMIT THE ROWS IN A.
C                   SIZE = N+1.
C NVA     JA    - COLUMN NUMBERS CORRESPONDING TO THE ELEMENTS OF A.
C                   SIZE = SIZE OF A.
C FVA     B     - RIGHT-HAND SIDE B;  B AND Z CAN THE SAME ARRAY.
C                   SIZE = N.
C FRA     Z     - SOLUTION X;  B AND Z CAN BE THE SAME ARRAY.
C                   SIZE = N.
C
C         THE ROWS AND COLUMNS OF THE ORIGINAL MATRIX M CAN BE
C    REORDERED (E.G., TO REDUCE FILLIN OR ENSURE NUMERICAL STABILITY)
C    BEFORE CALLING THE DRIVER.  IF NO REORDERING IS DONE, THEN SET
C    R(I) = C(I) = IC(I) = I  FOR I=1,...,N.  THE SOLUTION Z IS RETURNED
C    IN THE ORIGINAL ORDER.
C
C NVA     R     - ORDERING OF THE ROWS OF M.
C                   SIZE = N.
C NVA     C     - ORDERING OF THE COLUMNS OF M.
C                   SIZE = N.
C NVA     IC    - INVERSE OF THE ORDERING OF THE COLUMNS OF M;  I.E.,
C                   IC(C(I)) = I  FOR I=1,...,N.
C                   SIZE = N.
C
C         THE SOLUTION OF THE SYSTEM OF LINEAR EQUATIONS IS DIVIDED INTO
C    THREE STAGES --
C      NSF -- THE MATRIX M IS PROCESSED SYMBOLICALLY TO DETERMINE WHERE
C              FILLIN WILL OCCUR DURING THE NUMERIC FACTORIZATION.
C      NNF -- THE MATRIX M IS FACTORED NUMERICALLY INTO THE PRODUCT LDU
C              OF A UNIT LOWER TRIANGULAR MATRIX L, A DIAGONAL MATRIX D,
C              AND A UNIT UPPER TRIANGULAR MATRIX U, AND THE SYSTEM
C              MX = B  IS SOLVED.
C      NNS -- THE LINEAR SYSTEM  MX = B  IS SOLVED USING THE LDU
C  OR          FACTORIZATION FROM NNF.
C      NNT -- THE TRANSPOSED LINEAR SYSTEM  MT X = B  IS SOLVED USING
C              THE LDU FACTORIZATION FROM NNF.
C    FOR SEVERAL SYSTEMS WHOSE COEFFICIENT MATRICES HAVE THE SAME
C    NONZERO STRUCTURE, NSF NEED BE DONE ONLY ONCE (FOR THE FIRST
C    SYSTEM);  THEN NNF IS DONE ONCE FOR EACH ADDITIONAL SYSTEM.  FOR
C    SEVERAL SYSTEMS WITH THE SAME COEFFICIENT MATRIX, NSF AND NNF NEED
C    BE DONE ONLY ONCE (FOR THE FIRST SYSTEM);  THEN NNS OR NNT IS DONE
C    ONCE FOR EACH ADDITIONAL RIGHT-HAND SIDE.
C
C NV      PATH  - PATH SPECIFICATION;  VALUES AND THEIR MEANINGS ARE --
C                   1  PERFORM NSF AND NNF.
C                   2  PERFORM NNF ONLY  (NSF IS ASSUMED TO HAVE BEEN
C                       DONE IN A MANNER COMPATIBLE WITH THE STORAGE
C                       ALLOCATION USED IN THE DRIVER).
C                   3  PERFORM NNS ONLY  (NSF AND NNF ARE ASSUMED TO
C                       HAVE BEEN DONE IN A MANNER COMPATIBLE WITH THE
C                       STORAGE ALLOCATION USED IN THE DRIVER).
C                   4  PERFORM NNT ONLY  (NSF AND NNF ARE ASSUMED TO
C                       HAVE BEEN DONE IN A MANNER COMPATIBLE WITH THE
C                       STORAGE ALLOCATION USED IN THE DRIVER).
C                   5  PERFORM NSF ONLY.
C
C         VARIOUS ERRORS ARE DETECTED BY THE DRIVER AND THE INDIVIDUAL
C    SUBROUTINES.
C
C NR      FLAG  - ERROR FLAG;  VALUES AND THEIR MEANINGS ARE --
C                     0     NO ERRORS DETECTED
C                     N+K   NULL ROW IN A  --  ROW = K
C                    2N+K   DUPLICATE ENTRY IN A  --  ROW = K
C                    3N+K   INSUFFICIENT STORAGE IN NSF  --  ROW = K
C                    4N+1   INSUFFICIENT STORAGE IN NNF
C                    5N+K   NULL PIVOT  --  ROW = K
C                    6N+K   INSUFFICIENT STORAGE IN NSF  --  ROW = K
C                    7N+1   INSUFFICIENT STORAGE IN NNF
C                    8N+K   ZERO PIVOT  --  ROW = K
C                   10N+1   INSUFFICIENT STORAGE IN NDRV
C                   11N+1   ILLEGAL PATH SPECIFICATION
C
C         WORKING STORAGE IS NEEDED FOR THE FACTORED FORM OF THE MATRIX
C    M PLUS VARIOUS TEMPORARY VECTORS.  THE ARRAYS ISP AND RSP SHOULD BE
C    EQUIVALENCED;  INTEGER STORAGE IS ALLOCATED FROM THE BEGINNING OF
C    ISP AND REAL STORAGE FROM THE END OF RSP.
C
C NV      NSP   - DECLARED DIMENSION OF RSP;  NSP GENERALLY MUST
C                   BE LARGER THAN  5N+3 + 2K  (WHERE  K = (NUMBER OF
C                   NONZERO ENTRIES IN M)).
C NVIRA   ISP   - INTEGER WORKING STORAGE DIVIDED UP INTO VARIOUS ARRAYS
C                   NEEDED BY THE SUBROUTINES;  ISP AND RSP SHOULD BE
C                   EQUIVALENCED.
C                   SIZE = LRATIO*NSP
C FVIRA   RSP   - REAL WORKING STORAGE DIVIDED UP INTO VARIOUS ARRAYS
C                   NEEDED BY THE SUBROUTINES;  ISP AND RSP SHOULD BE
C                   EQUIVALENCED.
C                   SIZE = NSP.
C NR      ESP   - IF SUFFICIENT STORAGE WAS AVAILABLE TO PERFORM THE
C                   SYMBOLIC FACTORIZATION (NSF), THEN ESP IS SET TO THE
C                   AMOUNT OF EXCESS STORAGE PROVIDED (NEGATIVE IF
C                   INSUFFICIENT STORAGE WAS AVAILABLE TO PERFORM THE
C                   NUMERIC FACTORIZATION (NNF)).
C
C
C  CONVERSION TO DOUBLE PRECISION
C
C    TO CONVERT THESE ROUTINES FOR DOUBLE PRECISION ARRAYS, SIMPLY USE
C    THE DOUBLE PRECISION DECLARATIONS IN PLACE OF THE REAL DECLARATIONS
C    IN EACH SUBPROGRAM;  IN ADDITION, THE DATA VALUE OF THE INTEGER
C    VARIABLE LRATIO MUST BE SET AS INDICATED IN SUBROUTINE NDRV
C
        INTEGER  R(1), C(1), IC(1),  IA(1), JA(1),  ISP(1), ESP,
     *     PATH, FLAG,  Q, IM, D, U, ROW, TMP,  UMAX
C       REAL  A(1),  B(1),  Z(1),  RSP(1)
        DOUBLE PRECISION  A(1),  B(1),  Z(1),  RSP(1)
C
C  SET LRATIO EQUAL TO THE RATIO BETWEEN THE LENGTH OF FLOATING POINT
C  AND INTEGER ARRAY DATA;  E. G., LRATIO = 1 FOR (REAL, INTEGER),
C  LRATIO = 2 FOR (DOUBLE PRECISION, INTEGER)
C
        DATA LRATIO/2/
C
        IF (PATH.LT.1 .OR. 5.LT.PATH)  GO TO 111
C  ******  INITIALIZE AND DIVIDE UP TEMPORARY STORAGE  *****************
        IL = 1
        IU = IL + N+1
        JL = IU + N+1
C
C  ******  CALL NSF IF FLAG IS SET  ************************************
        IF ((PATH-1) * (PATH-5) .NE. 0)  GO TO 2
          MAX = (LRATIO*NSP + 1 - JL) - (N+1) - N
          JLMAX = MAX/2
          Q     = JL  + JLMAX
          IM    = Q   + (N+1)
          JUTMP = IM  +   N
          JUMAX = LRATIO*NSP + 1 - JUTMP
          ESP = MAX/LRATIO
          IF (JLMAX.LE.0 .OR. JUMAX.LE.0)  GO TO 110
          CALL  NSF8
     *       (N,  R, IC,  IA, JA,
     *        ISP(IL), ISP(JL), JLMAX,  ISP(IU), ISP(JUTMP), JUMAX,
     *        ISP(Q),  ISP(IM),  FLAG)
          IF (FLAG.NE.0)  GO TO 100
C  ******  MOVE JU NEXT TO JL  *****************************************
          JLMAX = ISP(IL+N)-1
          JU    = JL + JLMAX
          JUMAX = ISP(IU+N)-1
          IF (JUMAX.LE.0)  GO TO 2
          DO 1 J=1,JUMAX
   1        ISP(JU+J-1) = ISP(JUTMP+J-1)
C
C  ******  CALL REMAINING SUBROUTINES  *********************************
   2    JLMAX = ISP(IL+N)-1
        JU    = JL  + JLMAX
        JUMAX = ISP(IU+N)-1
        L     = (JU + JUMAX - 2 + LRATIO)  /  LRATIO    +    1
        LMAX  = JLMAX
        D     = L   + LMAX
        U     = D   + N
        ROW   = NSP + 1 - N
        TMP   = ROW - N
        UMAX  = TMP - U
        ESP = UMAX - JUMAX
C
        IF ((PATH-1) * (PATH-2) .NE. 0)  GO TO 3
          IF (UMAX.LE.0)  GO TO 110
          CALL  NNF8
     *       (N,  R, C, IC,  IA, JA, A,  Z,  B,
     *        ISP(IL), ISP(JL), RSP(L), LMAX,   RSP(D),
     *           ISP(IU), ISP(JU), RSP(U), UMAX,
     *        RSP(ROW),  RSP(TMP),  FLAG)
          IF (FLAG.NE.0)  GO TO 100
          RETURN
C
   3    IF ((PATH-3) .NE. 0)  GO TO 4
          CALL  NNS8
     *       (N,  R, C,
     *        ISP(IL), ISP(JL), RSP(L),  RSP(D),
     *          ISP(IU), ISP(JU), RSP(U),
     *        Z,  B,  RSP(TMP))
C
   4    IF ((PATH-4) .NE. 0)  GO TO 5
          CALL  NNT8
     *       (N,  R, C,
     *        ISP(IL), ISP(JL), RSP(L),  RSP(D),
     *          ISP(IU), ISP(JU), RSP(U),
     *        Z,  B,  RSP(TMP))
   5    RETURN
C
C ** ERROR:  ERROR DETECTED IN NSF, NNF, NNS, OR NNT
 100    RETURN
C ** ERROR:  INSUFFICIENT STORAGE
 110    FLAG = 10*N + 1
        RETURN
C ** ERROR:  ILLEGAL PATH SPECIFICATION
 111    FLAG = 11*N + 1
        RETURN
        END
C
C       ----------------------------------------------------------------
C
C               YALE SPARSE MATRIX PACKAGE - NONSYMMETRIC CODES
C                    SOLVING THE SYSTEM OF EQUATIONS MX = B
C                        (UNCOMPRESSED POINTER STORAGE)
C
C    I.   CALLING SEQUENCES
C         THE COEFFICIENT MATRIX CAN BE PROCESSED BY AN ORDERING ROUTINE
C    (E.G., TO REDUCE FILLIN OR ENSURE NUMERICAL STABILITY) BEFORE USING
C    THE REMAINING SUBROUTINES.  IF NO REORDERING IS DONE, THEN SET
C    R(I) = C(I) = IC(I) = I  FOR I=1,...,N.  THE CALLING SEQUENCE IS --
C        (      (MATRIX ORDERING))
C         NSF   (SYMBOLIC FACTORIZATION TO DETERMINE WHERE FILLIN WILL
C                 OCCUR DURING NUMERIC FACTORIZATION)
C         NNF   (NUMERIC FACTORIZATION INTO PRODUCT LDU OF UNIT LOWER
C                 TRIANGULAR MATRIX L, DIAGONAL MATRIX D, AND UNIT UPPER
C                 TRIANGULAR MATRIX U, AND SOLUTION OF LINEAR SYSTEM)
C         NNS   (SOLUTION OF LINEAR SYSTEM FOR ADDITIONAL RIGHT-HAND
C     OR          SIDE USING LDU FACTORIZATION FROM NNF)
C         NNT   (SOLUTION OF TRANSPOSED LINEAR SYSTEM FOR ADDITIONAL
C                 RIGHT-HAND SIDE USING LDU FACTORIZATION FROM NNF)
C
C    II.  STORAGE OF SPARSE MATRICES
C         THE NONZERO ENTRIES OF THE COEFFICIENT MATRIX M ARE STORED
C    ROW-BY-ROW IN THE ARRAY A.  TO IDENTIFY THE INDIVIDUAL NONZERO
C    ENTRIES IN EACH ROW, WE NEED TO KNOW IN WHICH COLUMN EACH ENTRY
C    LIES.  THE COLUMN INDICES WHICH CORRESPOND TO THE NONZERO ENTRIES
C    OF M ARE STORED IN THE ARRAY JA;  I.E., IF  A(K) = M(I,J),  THEN
C    JA(K) = J.  IN ADDITION, WE NEED TO KNOW WHERE EACH ROW STARTS AND
C    HOW LONG IT IS.  THE INDEX POSITIONS IN JA AND A WHERE THE ROWS OF
C    M BEGIN ARE STORED IN THE ARRAY IA;  I.E., IF M(I,J) IS THE FIRST
C    NONZERO ENTRY (STORED) IN THE I-TH ROW AND A(K) = M(I,J),  THEN
C    IA(I) = K.  MOREOVER, THE INDEX IN JA AND A OF THE FIRST LOCATION
C    FOLLOWING THE LAST ELEMENT IN THE LAST ROW IS STORED IN IA(N+1).
C    THUS, THE NUMBER OF ENTRIES IN THE I-TH ROW IS GIVEN BY
C    IA(I+1) - IA(I),  THE NONZERO ENTRIES OF THE I-TH ROW ARE STORED
C    CONSECUTIVELY IN
C            A(IA(I)),  A(IA(I)+1),  ..., A(IA(I+1)-1),
C    AND THE CORRESPONDING COLUMN INDICES ARE STORED CONSECUTIVELY IN
C            JA(IA(I)), JA(IA(I)+1), ..., JA(IA(I+1)-1).
C    FOR EXAMPLE, THE 5 BY 5 MATRIX
C                ( 1. 0. 2. 0. 0.)
C                ( 0. 3. 0. 0. 0.)
C            M = ( 0. 4. 5. 6. 0.)
C                ( 0. 0. 0. 7. 0.)
C                ( 0. 0. 0. 8. 9.)
C    WOULD BE STORED AS
C                 1  2  3  4  5  6  7  8  9
C            ---+--------------------------
C            IA   1  3  4  7  8 10
C            JA   1  3  2  2  3  4  4  4  5
C             A   1. 2. 3. 4. 5. 6. 7. 8. 9.         .
C
C         THE STRICT TRIANGULAR PORTIONS OF THE MATRICES L AND U ARE
C    STORED IN THE SAME FASHION USING THE ARRAYS  IL, JL, L  AND
C    IU, JU, U  RESPECTIVELY.  THE DIAGONAL ENTRIES OF L AND U ARE
C    ASSUMED TO BE EQUAL TO ONE AND ARE NOT STORED.  THE ARRAY D
C    CONTAINS THE RECIPROCALS OF THE DIAGONAL ENTRIES OF THE MATRIX D.
C
C    III. ADDITIONAL STORAGE SAVINGS
C         IN NSF, R AND IC CAN BE THE SAME ARRAY IN THE CALLING
C    SEQUENCE IF NO REORDERING OF THE COEFFICIENT MATRIX HAS BEEN DONE.
C         IN NNF, R, C AND IC CAN ALL BE THE SAME ARRAY IF NO REORDERING
C    HAS BEEN DONE.  IF ONLY THE ROWS HAVE BEEN REORDERED, THEN C AND IC
C    CAN BE THE SAME ARRAY.  IF THE ROW AND COLUMN ORDERINGS ARE THE
C    SAME, THEN R AND C CAN BE THE SAME ARRAY.  Z AND ROW CAN BE THE
C    SAME ARRAY.
C         IN NNS OR NNT, R AND C CAN BE THE SAME ARRAY IF NO REORDERING
C    HAS BEEN DONE OR IF THE ROW AND COLUMN ORDERINGS ARE THE SAME.  Z
C    AND B CAN BE THE SAME ARRAY;  HOWEVER, THEN B WILL BE DESTROYED.
C
C    IV.  PARAMETERS
C         FOLLOWING IS A LIST OF PARAMETERS TO THE PROGRAMS.  NAMES ARE
C    UNIFORM AMONG THE VARIOUS SUBROUTINES.  CLASS ABBREVIATIONS ARE --
C       N - INTEGER VARIABLE
C       F - REAL VARIABLE
C       V - SUPPLIES A VALUE TO A SUBROUTINE
C       R - RETURNS A RESULT FROM A SUBROUTINE
C       I - USED INTERNALLY BY A SUBROUTINE
C       A - ARRAY
C
C CLASS   PARAMETER
C ------+----------
C FVA     A     - NONZERO ENTRIES OF THE COEFFICIENT MATRIX M, STORED
C                   BY ROWS.
C                   SIZE = NUMBER OF NONZERO ENTRIES IN M.
C FVA     B     - RIGHT-HAND SIDE B.
C                   SIZE = N.
C NVA     C     - ORDERING OF THE COLUMNS OF M.
C                   SIZE = N.
C FVRA    D     - RECIPROCALS OF THE DIAGONAL ENTRIES OF THE MATRIX D.
C                   SIZE = N.
C NR      FLAG  - ERROR FLAG;  VALUES AND THEIR MEANINGS ARE --
C                    0     NO ERRORS DETECTED
C                    N+K   NULL ROW IN A  --  ROW = K
C                   2N+K   DUPLICATE ENTRY IN A  --  ROW = K
C                   3N+K   INSUFFICIENT STORAGE FOR JL  --  ROW = K
C                   4N+1   INSUFFICIENT STORAGE FOR L
C                   5N+K   NULL PIVOT  --  ROW = K
C                   6N+K   INSUFFICIENT STORAGE FOR JU  --  ROW = K
C                   7N+1   INSUFFICIENT STORAGE FOR U
C                   8N+K   ZERO PIVOT  --  ROW = K
C NVA     IA    - POINTERS TO DELIMIT THE ROWS IN A.
C                   SIZE = N+1.
C NVA     IC    - INVERSE OF THE ORDERING OF THE COLUMNS OF M;  I.E.,
C                   IC(C(I) = I  FOR I=1,...N.
C                   SIZE = N.
C NVRA    IL    - POINTERS TO DELIMIT THE ROWS IN L.
C                   SIZE = N+1.
C NVRA    IU    - POINTERS TO DELIMIT THE ROWS IN U.
C                   SIZE = N+1.
C NVA     JA    - COLUMN NUMBERS CORRESPONDING TO THE ELEMENTS OF A.
C                   SIZE = SIZE OF A.
C NVRA    JL    - COLUMN NUMBERS CORRESPONDING TO THE ELEMENTS OF L.
C                   SIZE = JLMAX.
C NV      JLMAX - DECLARED DIMENSION OF JL;  JLMAX MUST BE LARGER THAN
C                   THE NUMBER OF NONZERO ENTRIES IN THE STRICT LOWER
C                   TRIANGLE OF M PLUS FILLIN  (IL(N+1)-1 AFTER NSF).
C NVRA    JU    - COLUMN NUMBERS CORRESPONDING TO THE ELEMENTS OF U.
C                   SIZE = JUMAX.
C NV      JUMAX - DECLARED DIMENSION OF JU;  JUMAX MUST BE LARGER THAN
C                   THE NUMBER OF NONZERO ENTRIES IN THE STRICT UPPER
C                   TRIANGLE OF M PLUS FILLIN  (IU(N+1)-1 AFTER NSF).
C FVRA    L     - NONZERO ENTRIES IN THE STRICT LOWER TRIANGULAR PORTION
C                   OF THE MATRIX L, STORED BY ROWS.
C                   SIZE = LMAX
C NV      LMAX  - DECLARED DIMENSION OF L;  LMAX MUST BE LARGER THAN
C                   THE NUMBER OF NONZERO ENTRIES IN THE STRICT LOWER
C                   TRIANGLE OF M PLUS FILLIN  (IL(N+1)-1 AFTER NSF).
C NV      N     - NUMBER OF VARIABLES/EQUATIONS.
C NVA     R     - ORDERING OF THE ROWS OF M.
C                   SIZE = N.
C FVRA    U     - NONZERO ENTRIES IN THE STRICT UPPER TRIANGULAR PORTION
C                   OF THE MATRIX U, STORED BY ROWS.
C                   SIZE = UMAX.
C NV      UMAX  - DECLARED DIMENSION OF U;  UMAX MUST BE LARGER THAN
C                   THE NUMBER OF NONZERO ENTRIES IN THE STRICT UPPER
C                   TRIANGLE OF M PLUS FILLIN  (IU(N+1)-1 AFTER NSF).
C FRA     Z     - SOLUTION X.
C                   SIZE = N.
C
C
C       ----------------------------------------------------------------
C
C*** SUBROUTINE NSF
C*** SYMBOLIC LDU-FACTORIZATION OF A NONSYMMETRIC SPARSE MATRIX
C      (UNCOMPRESSED POINTER STORAGE)
C
        SUBROUTINE  NSF8
     *     (N, R,IC, IA,JA, IL,JL,JLMAX, IU,JU,JUMAX, Q, IM, FLAG)
C
C       INPUT VARIABLES:   N, R,IC, IA,JA, JLMAX, JUMAX.
C       OUTPUT VARIABLES:  IL,JL, IU,JU, FLAG.
C
C       PARAMETERS USED INTERNALLY:
C NIA     Q     - SUPPOSE M' IS THE RESULT OF REORDERING M;  IF
C                   PROCESSING OF THE KTH ROW OF M' (HENCE THE KTH ROWS
C                   OF L AND U) IS BEING DONE, THEN Q(J) IS INITIALLY
C                   NONZERO IF M'(K,J) IS NONZERO;  SINCE VALUES NEED
C                   NOT BE STORED, EACH ENTRY POINTS TO THE NEXT
C                   NONZERO;  FOR EXAMPLE, IF  N=9  AND THE 5TH ROW OF
C                   M' IS
C                           0 X X 0 X 0 0 X 0,
C                   THEN Q WILL INITIALLY BE
C                           A 3 5 A 8 A A 10 A 2        (A - ARBITRARY);
C                   Q(N+1) POINTS TO THE FIRST NONZERO IN THE ROW AND
C                   THE LAST NONZERO POINTS TO  N+1;  AS THE ALGORITHM
C                   PROCEEDS, OTHER ELEMENTS OF Q ARE INSERTED IN THE
C                   LIST BECAUSE OF FILLIN.
C                   SIZE = N+1.
C NIA     IM    - AT EACH STEP IN THE FACTORIZATION, IM(I) IS THE LAST
C                   ELEMENT IN THE ITH ROW OF U WHICH NEEDS TO BE
C                   CONSIDERED IN COMPUTING FILLIN.
C                   SIZE = N.
C
C  INTERNAL VARIABLES--
C    JLPTR - POINTS TO THE LAST POSITION USED IN  JL.
C    JUPTR - POINTS TO THE LAST POSITION USED IN  JU.
C
        INTEGER  R(1), IC(1),  IA(1), JA(1),  IL(1), JL(1),
     *     IU(1), JU(1),  Q(1),  IM(1),  FLAG,  QM, VJ
C
C  ******  INITIALIZE POINTERS  ****************************************
        JLPTR = 0
        IL(1) = 1
        JUPTR = 0
        IU(1) = 1
C
C  ******  FOR EACH ROW OF L AND U  ************************************
        DO 10 K=1,N
C  ******  SET Q TO THE REORDERED ROW OF A  ****************************
          Q(N+1) = N+1
          JMIN = IA(R(K))
          JMAX = IA(R(K)+1) - 1
          IF (JMIN.GT.JMAX)  GO TO 101
          DO 2 J=JMIN,JMAX
            VJ = IC(JA(J))
            QM = N+1
   1        M = QM
            QM = Q(M)
            IF (QM.LT.VJ)  GO TO 1
            IF (QM.EQ.VJ)  GO TO 102
              Q(M) = VJ
              Q(VJ) = QM
   2        CONTINUE
C
C  ******  FOR EACH ENTRY IN THE LOWER TRIANGLE  ***********************
          I = N+1
   3      I = Q(I)
          IF (I.GE.K)  GO TO 7
C  ******  L(K,I) WILL BE NONZERO, SO ADD IT TO JL  ********************
            JLPTR = JLPTR+1
            IF (JLPTR.GT.JLMAX)  GO TO 103
            JL(JLPTR) = I
            QM = I
C  ******  INSPECT ITH ROW FOR FILLIN, ADJUST IM IF POSSIBLE  **********
            JMIN = IU(I)
            JMAX = IM(I)
            IF (JMIN.GT.JMAX)  GO TO 6
            DO 5 J=JMIN,JMAX
              VJ = JU(J)
              IF (VJ.EQ.K)  IM(I) = J
   4          M = QM
              QM = Q(M)
              IF (QM.LT.VJ)  GO TO 4
              IF (QM.EQ.VJ)  GO TO 5
                Q(M) = VJ
                Q(VJ) = QM
                QM = VJ
   5          CONTINUE
   6        GO TO 3
C
C  ******  CHECK FOR NULL PIVOT  ***************************************
   7      IF (I.NE.K)  GO TO 105
C  ******  REMAINING ELEMENTS OF Q DEFINE STRUCTURE OF U(K, )  *********
   8      I = Q(I)
          IF (I.GT.N)  GO TO 9
            JUPTR = JUPTR+1
            IF (JUPTR.GT.JUMAX)  GO TO 106
            JU(JUPTR) = I
            GO TO 8
C  ******  GET READY FOR NEXT ROW  *************************************
   9      IM(K) = JUPTR
          IL(K+1) = JLPTR+1
  10      IU(K+1) = JUPTR+1
C
        FLAG = 0
        RETURN
C
C ** ERROR:  NULL ROW IN A
 101    FLAG = N + R(K)
        RETURN
C ** ERROR:  DUPLICATE ENTRY IN A
 102    FLAG = 2*N + R(K)
        RETURN
C ** ERROR:  INSUFFICIENT STORAGE FOR JL
 103    FLAG = 3*N + K
        RETURN
C ** ERROR:  NULL PIVOT
 105    FLAG = 5*N + K
        RETURN
C ** ERROR:  INSUFFICIENT STORAGE FOR JU
 106    FLAG = 6*N + K
        RETURN
        END
C
C       ----------------------------------------------------------------
C
C*** SUBROUTINE NNF
C*** NUMERIC LDU-FACTORIZATION OF SPARSE NONSYMMETRIC MATRIX AND
C      SOLUTION OF SYSTEM OF LINEAR EQUATIONS (UNCOMPRESSED POINTER
C      STORAGE)
C
        SUBROUTINE  NNF8
     *     (N, R,C,IC, IA,JA,A, Z, B, IL,JL,L,LMAX, D, IU,JU,U,UMAX,
     *      ROW, TMP, FLAG)
C
C       INPUT VARIABLES:   N, R,C,IC, IA,JA,A, B, IL,JL,LMAX, IU,JU,UMAX
C       OUTPUT VARIABLES:  Z, L,D,U, FLAG
C
C       PARAMETERS USED INTERNALLY:
C FIA     ROW   - HOLDS INTERMEDIATE VALUES IN CALCULATION OF L, D, U.
C                   SIZE = N.
C FIA     TMP   - HOLDS NEW RIGHT-HAND SIDE B' FOR SOLUTION OF THE
C                   EQUATION  UX = B'.
C                   SIZE = N.
C
        INTEGER  R(1), C(1), IC(1),  IA(1), JA(1),
     *     IL(1), JL(1), LMAX,  IU(1), JU(1), UMAX,  FLAG
C       REAL  A(1), Z(1), B(1),  L(1), D(1), U(1),
C    *     ROW(1), TMP(1),  LI, SUM, DK
        DOUBLE PRECISION  A(1), Z(1), B(1),  L(1), D(1), U(1),
     *     ROW(1), TMP(1),  LI, SUM, DK
C
C  ******  CHECK STORAGE  **********************************************
        IF (IL(N+1)-1 .GT. LMAX)  GO TO 104
        IF (IU(N+1)-1 .GT. UMAX)  GO TO 107
C
C  ******  FOR EACH ROW  ***********************************************
        DO 10 K=1,N
C  ******  SET THE INITIAL STRUCTURE OF ROW  ***************************
          JMIN = IL(K)
          JMAX = IL(K+1) - 1
          IF (JMIN.GT.JMAX)  GO TO 2
C  ******  IF L(K,M) .NE. 0, ROW(M)=0  *********************************
          DO 1 J=JMIN,JMAX
   1        ROW(JL(J)) = 0
   2      ROW(K) = 0
          JMIN = IU(K)
          JMAX = IU(K+1) - 1
          IF (JMIN.GT.JMAX)  GO TO 4
C  ******  IF U(K,M) .NE. 0, ROW(M)=0  *********************************
          DO 3 J=JMIN,JMAX
   3        ROW(JU(J)) = 0
   4      JMIN = IA(R(K))
          JMAX = IA(R(K)+1) - 1
C  ******  SET ROW TO KTH ROW OF REORDERED A  **************************
          DO 5 J=JMIN,JMAX
   5        ROW(IC(JA(J))) = A(J)
C  ******  INITIALIZE SUM  *********************************************
          SUM = B(R(K))
C
C  ******  ASSIGN THE KTH ROW OF L AND ADJUST ROW, SUM  ****************
          IMIN = IL(K)
          IMAX = IL(K+1) - 1
          IF (IMIN.GT.IMAX)  GO TO 8
          DO 7 I=IMIN,IMAX
            LI = - ROW(JL(I))
C  ******  IF L IS NOT REQUIRED, THEN COMMENT OUT THE FOLLOWING LINE  **
            L(I) = - LI
            SUM = SUM + LI * TMP(JL(I))
            JMIN = IU(JL(I))
            JMAX = IU(JL(I)+1) - 1
            IF (JMIN.GT.JMAX)  GO TO 7
            DO 6 J=JMIN,JMAX
   6          ROW(JU(J)) = ROW(JU(J)) + LI * U(J)
   7        CONTINUE
C
C  ******  ASSIGN DIAGONAL D AND KTH ROW OF U, SET TMP(K)  *************
   8      IF (ROW(K).EQ.0)  GO TO 108
          DK = 1 / ROW(K)
          D(K) = DK
          TMP(K) = SUM * DK
          JMIN = IU(K)
          JMAX = IU(K+1) - 1
          IF (JMIN.GT.JMAX)  GO TO 10
          DO 9 J=JMIN,JMAX
   9        U(J) = ROW(JU(J)) * DK
  10      CONTINUE
C
C  ******  SOLVE  UX = TMP  BY BACK SUBSTITUTION  **********************
        K = N
        DO 13 I=1,N
          SUM = TMP(K)
          JMIN = IU(K)
          JMAX = IU(K+1) - 1
          IF (JMIN.GT.JMAX)  GO TO 12
          DO 11 J=JMIN,JMAX
  11        SUM = SUM - U(J) * TMP(JU(J))
  12      TMP(K) = SUM
          Z(C(K)) = SUM
  13      K = K-1
C
        FLAG = 0
        RETURN
C
C ** ERROR:  INSUFFICIENT STORAGE FOR L
 104    FLAG = 4*N + 1
        RETURN
C ** ERROR:  INSUFFICIENT STORAGE FOR U
 107    FLAG = 7*N + 1
        RETURN
C ** ERROR:  ZERO PIVOT
 108    FLAG = 8*N + K
        RETURN
        END
C
C       ----------------------------------------------------------------
C
C*** SUBROUTINE NNS
C*** NUMERIC SOLUTION OF A SPARSE NONSYMMETRIC SYSTEM OF LINEAR
C      EQUATIONS GIVEN LDU-FACTORIZATION (UNCOMPRESSED POINTER STORAGE)
C
        SUBROUTINE  NNS8
     *     (N, R,C, IL,JL,L, D, IU,JU,U, Z, B, TMP)
C
C       INPUT VARIABLES:   N, R,C, IL,JL,L, D, IU,JU,U, B
C       OUTPUT VARIABLES:  Z
C
C       PARAMETERS USED INTERNALLY:
C FIA     TMP   - HOLDS NEW RIGHT-HAND SIDE B' FOR SOLUTION OF THE
C                   EQUATION UX = B'.
C                   SIZE = N.
C
        INTEGER  R(1), C(1),  IL(1), JL(1),  IU(1), JU(1)
C       REAL  L(1), D(1), U(1),  Z(1), B(1),  TMP(1), SUM
        DOUBLE PRECISION  L(1), D(1), U(1),  Z(1), B(1),  TMP(1), SUM
C
C  ******  SOLVE LDY = B  BY FORWARD SUBSTITUTION  *********************
        DO 2 K=1,N
          SUM = B(R(K))
          JMIN = IL(K)
          JMAX = IL(K+1) - 1
          IF (JMIN.GT.JMAX)  GO TO 2
          DO 1 J=JMIN,JMAX
   1        SUM = SUM - L(J) * TMP(JL(J))
   2      TMP(K) = SUM * D(K)
C
C  ******  SOLVE  UX = Y  BY BACK SUBSTITUTION  ************************
        K = N
        DO 5 I=1,N
          SUM = TMP(K)
          JMIN = IU(K)
          JMAX = IU(K+1) - 1
          IF (JMIN.GT.JMAX)  GO TO 4
          DO 3 J=JMIN,JMAX
   3        SUM = SUM - U(J) * TMP(JU(J))
   4      TMP(K) = SUM
          Z(C(K)) = SUM
   5      K = K-1
        RETURN
        END
C       ----------------------------------------------------------------
C
C*** SUBROUTINE NNT
C*** NUMERIC SOLUTION OF THE TRANSPOSE OF A SPARSE NONSYMMETRIC SYSTEM
C      OF LINEAR EQUATIONS GIVEN LDU-FACTORIZATION (UNCOMPRESSED POINTER
C      STORAGE)
C
        SUBROUTINE  NNT8
     *     (N, R,C, IL,JL,L, D, IU,JU,U, Z, B, TMP)
C
C       INPUT VARIABLES:   N, R,C, IL,JL,L, D, IU,JU,U, B
C       OUTPUT VARIABLES:  Z
C
C       PARAMETERS USED INTERNALLY:
C FIA     TMP   - HOLDS NEW RIGHT-HAND SIDE B' FOR SOLUTION OF THE
C                   EQUATION LX = B'.
C                   SIZE = N.
C
        INTEGER  R(1), C(1),  IL(1), JL(1),  IU(1), JU(1)
C       REAL  L(1), D(1), U(1),  Z(1), B(1),  TMP(1), TMPK
        DOUBLE PRECISION  L(1), D(1), U(1),  Z(1), B(1),  TMP(1), TMPK
C
C  ******  SOLVE  UT Y = B  BY FORWARD SUBSTITUTION  *******************
        DO 1 K=1,N
   1      TMP(K) = B(C(K))
        DO 3 K=1,N
          TMPK = - TMP(K)
          JMIN = IU(K)
          JMAX = IU(K+1) - 1
          IF (JMIN.GT.JMAX)  GO TO 3
          DO 2 J=JMIN,JMAX
   2        TMP(JU(J)) = TMP(JU(J)) + U(J) * TMPK
   3      CONTINUE
C
C  ******  SOLVE  D LT X = Y  BY BACK SUBSTITUTION  ********************
        K = N
        DO 6 I=1,N
          TMPK = - TMP(K) * D(K)
          JMIN = IL(K)
          JMAX = IL(K+1) - 1
          IF (JMIN.GT.JMAX)  GO TO 5
          DO 4 J=JMIN,JMAX
   4        TMP(JL(J)) = TMP(JL(J)) + L(J) * TMPK
   5      Z(R(K)) = - TMPK
   6      K = K-1
        RETURN
        END
